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Abstract 

Stochastic  particle  methods  (SPMs)  for  the  Boltzmann  equation,  such  as  the  Direct 
Simulation  Monte  Carlo  (DSMC)  technique,  have  gained  popularity  for  the  prediction 
of  flows  in  which  the  assumptions  behind  the  continuum  equations  of  fluid  mechanics 
break  down;  however,  there  are  still  a  number  of  issues  that  make  SPMs  computationally 
challenging  for  practical  use.  In  traditional  SPMs,  simulated  particles  may  possess  only 
a  single  velocity  vector,  even  though  they  may  represent  an  extremely  large  collection 
of  actual  particles.  This  limits  the  method  to  converge  only  in  law  to  the  Boltzmann 
solution.  This  document  details  the  development  of  new  SPMs  that  allow  the  velocity 
of  each  simulated  particle  to  be  distributed.  This  approach  has  been  termed  Distributional 
Monte  Carlo  (DMC). 

A  technique  is  described  which  applies  kernel  density  estimation  to  Nanbu’s  DSMC 
algorithm.  It  is  then  proven  that  the  method  converges  not  just  in  law,  but  also  in  solution 
for  L°°  (M3J  solutions  of  the  space  homogeneous  Boltzmann  equation.  This  provides  for 
direct  evaluation  of  the  velocity  density  function.  The  derivation  of  a  general  Distributional 
Monte  Carlo  method  is  given  which  treats  collision  interactions  between  simulated  particles 
as  a  relaxation  problem.  The  framework  is  proven  to  converge  in  law  to  the  solution  of  the 
space  homogeneous  Boltzmann  equation,  as  well  as  in  solution  for  L°°  (M3)  solutions.  An 
approach  based  on  the  BGK  simplification  is  presented  which  computes  collision  outcomes 
deterministically. 

Each  technique  is  applied  to  the  well-studied  Bobylev-Krook-Wu  solution  as  a 
numerical  test  case.  Accuracy  and  variance  of  the  solutions  are  examined  as  functions 
of  various  simulation  parameters.  Significantly  improved  accuracy  and  reduced  variance 
are  observed  in  the  normalized  moments  for  the  Distributional  Monte  Carlo  technique 
employing  discrete  BGK  collision  modeling. 
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DISTRIBUTIONAL  MONTE  CARLO  METHODS  FOR 


THE  BOLTZMANN  EQUATION 

I.  Introduction 

The  study  of  the  thermodynamic  properties  of  gas  flows  is  critical  in  many  fields 
of  engineering  and  the  sciences.  The  field  of  computational  fluid  mechanics  has  in 
recent  years  made  ever  increasing  strides  in  the  analysis  of  fluids  governed  by  the 
continuum  equations  of  fluid  mechanics  (e.g.  the  Navier-Stokes  and  Euler  equations). 
Numerical  methods  for  solving  these  equation  sets  have  found  their  way  into  many  practical 
engineering  tools  and  practices,  and  many  recent  developments  in  applied  aerodynamics 
can  be  attributed  to  their  use. 

Unfortunately,  the  continuum  equation  sets  cannot  provide  a  complete  description 
of  the  physical  phenomena  taking  place  within  a  fluid  under  all  circumstances.  Notable 
departures  from  the  predicted  solutions  of  these  equations  occur  whenever  the  “continuum 
hypothesis”  is  violated.  Examples  of  such  violations  arise  in  rarefied  gas  dynamics, 
hypersonic  flows,  and  micro-scale  flows.  Additionally,  there  is  an  “equilibrium  hypothesis” 
inherent  to  the  so-called  “continuum”  equation  sets,  as  it  can  be  shown  that  these  equation 
sets  may  be  derived  by  assuming  various  forms  of  a  perturbation  from  local  thermodynamic 
equilibrium  [32]. 

To  obtain  accurate  solutions  in  such  cases,  models  based  on  kinetic  theory  are  used. 
Kinetic  theory  attempts  to  relate  the  molecular  interactions  occurring  at  the  microscopic 
level  to  macroscopic  fluid  properties  such  as  pressure,  temperature,  viscosity,  etc.  The  most 
common  governing  equation  employed  is  the  Boltzmann  equation,  which  is  an  integro- 
differential  equation  which  describes  the  evolution  of  the  probability  density  function 
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(PDF)  of  molecular  velocities  throughout  the  gas.  This  function  is  commonly  referred 
to  in  the  literature  as  the  velocity  distribution  function,  but  to  be  more  mathematically 
precise  this  document  will  utilize  the  term  velocity  density  function  when  referring  to  the 
Boltzmann  solution. 

Unlike  the  continuum  equation  sets,  few  practical  solution  methods  for  the  Boltzmann 
equation  have  found  their  way  into  practical  engineering  settings.  This  is  not  for  a  lack 
of  solution  algorithms,  but  rather  due  to  the  complexity  and  computational  demands  of 
such  algorithms.  By  far,  the  Direct  Simulation  Monte  Carlo  (DSMC)  method,  originally 
developed  by  Bird  [14],  has  become  the  dominant  approach  for  examining  flows  governed 
by  the  Boltzmann  equation  and  has  gained  general  acceptance  for  practical  applications  in 
rarefied  gas  dynamics. 

The  DSMC  technique  approximates  the  physics  of  the  Boltzmann  equation  using  a 
stochastic  simulation  of  the  interactions  of  a  fraction  of  the  molecules  in  the  gas.  In 
DSMC,  each  simulated  particle  possesses  a  single  velocity  vector  and  energy  state.  As 
only  a  fraction  of  the  number  of  particles  in  the  gas  can  be  simulated,  each  simulated 
particle  is  assumed  to  represent  W  =  actual  particles,  where  N  is  the  total  number  of 
actual  particles  in  the  gas  and  Np  is  the  number  of  simulated  particles.  In  practice  W  may  be 
quite  large  (106  or  more),  therefore  restricting  to  only  a  single  velocity  vector  per  simulated 
particle  is  equivalent  to  assuming  that  millions  of  actual  particles  all  share  the  exact  velocity 
vector.  This  assumption  is  non-physical  in  the  sense  that  kinetic  theory  tells  us  that  the 
probability  that  any  two  molecules  share  the  exact  velocity  is  zero.  Mathematically,  this 
representation  leads  to  what  is  known  as  a  point  measure  approximation  of  the  density 
function,  /,  as  illustrated  by 

1  np 

f(y)  =  —  Vd(v  — V/)  (1.1) 

lyp  i= 1 

where  v,-  is  the  velocity  vector  of  the  ith  simulated  particle,  and  6  is  the  Dirac  distribution. 
This  representation  restricts  the  method  to  converge  only  in  law.  Another  undesirable 
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feature  of  the  DSMC  method  stems  from  the  stochastic  nature  of  the  technique.  Namely, 
flow  solutions  are  subject  to  a  significant  amount  of  variance.  In  practice  this  variance  is 
reduced  by  creating  an  ensemble  averaged  solution,  combining  the  results  of  potentially 
thousands  of  simulations  (or  time  steps  in  steady-state  cases).  The  present  work  seeks  to 
address  these  limitations  through  the  development  of  Distributional  Monte  Carlo  methods. 
Although  the  Boltzmann  equation  and  DSMC  will  be  discussed  in  greater  detail  in  Chapters 
2  and  3,  a  brief  introduction  is  warranted  here  to  provide  a  more  comprehensive  view  of 
the  motivation  and  contribution  of  the  work  presented  herein. 

1.1  Kinetic  Theory  and  Rarefied  Gas  Dynamics 

Although  the  continuum  equations  are  capable  of  providing  fairly  accurate  flowfields 
under  relatively  benign  conditions,  two  specific  regimes  in  which  they  fail  are  rarefied  flows 
and  flows  containing  non-equilibrium  phenomena.  Applications  of  rarefied  gas  dynamics 
typically  involve  high-altitude  flight  and  microscale  flows  (e.g.  Micro  Electro-Mechanical 
devices(MEMs)).  In  the  former,  the  atmospheric  density  is  low  enough  that  the  large 
intermolecular  spacing  invalidates  the  continuum  hypothesis.  In  the  latter,  the  physical 
scale  of  interest  is  small  enough  that  the  flow  appears  rarefied  even  at  standard  densities. 

Non-equilibrium  phenomena  typically  result  from  the  propensity  of  the  constituent 
molecules  to  undergo  changes  in  their  internal  energy  state  during  a  collision  with  a 
material  surface  or  another  molecule,  as  well  as  their  ability  to  react  chemically  with 
other  molecules  upon  collision.  These  events  occur  at  some  finite  rate  in  the  fluid,  and 
not  every  encounter  results  in  such  changes.  If  these  events  occur  in  such  a  way  as  to  alter 
the  macroscopic  properties  of  the  fluid  over  timescales  which  are  longer  than  some  local 
flow  timescale,  they  may  be  regarded  as  non-equilibrium  phenomena.  These  effects  can 
be  responsible  for  causing  a  number  of  fluid  properties  that  might  normally  be  treated  as 
constants  to  vary  at  a  finite  rate.  Among  these  are  fluid  composition,  viscosity,  thermal 
conductivity,  specific  heats,  etc.  The  variable  nature  of  these  properties  in  regions  of  non- 
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equilibrium  changes  the  manner  in  which  energy  and  momentum  are  transferred  in  the 
fluid,  which  alters  the  macroscopic  thermodynamic  and  flow  variables  throughout  the  flow. 
Such  effects  are  commonly  observed  in  high  temperature  or  high  energy  flow  fields  such  as 
those  generated  in  hypersonic  flight. 

The  traditional  continuum  equations  of  fluid  mechanics  do  nothing  in  and  of 
themselves  to  address  these  effects.  In  fact,  these  equations  permit  only  small  departures 
from  equilibrium  in  the  translational  energy  mode,  but  must  be  augmented  with  additional 
equations  to  compensate  for  other  forms  of  non-equilibrium.  Integration  of  such  models 
with  the  continuum  equations  in  some  cases  is  not  well  understood. 

To  illustrate  the  need  for  a  noncontinuum  based  approach,  consider  a  blunt  body 
in  a  hypersonic  flow,  as  illustrated  in  Figure  1.1.  The  physics  involved  with  the 
associated  strong  shocks  invalidate  the  continuum  equations  through  the  generation  of  non¬ 
equilibrium  phenomena.  Close  to  the  body,  strong  gradients  exist  across  the  boundary  layer 
and  may  generate  non-equilibrium  phenomena  which  invalidate  the  continuum  equations. 
Further  downstream  the  flow  may  expand  around  the  body  to  the  point  at  which  the  density 
is  too  low  for  the  continuum  equations  to  hold.  Clearly  the  continuum  equations  are  not 
sufficient  for  dealing  with  such  flows. 


Figure  1.1:  Hypersonic  Flow  over  a  Blunt  Body 


4 


A  common  parameter  used  to  examine  the  validity  of  the  continuum  equations  is  the 


Knudsen  Number,  Kn,  defined  by 


(1.2) 


Here,  A  represents  the  mean  free  path  of  a  molecule  in  the  gas  and  L  is  a  characteristic 
length.  The  characteristic  length  of  importance  depends  upon  the  situation  under 
consideration.  For  example,  to  obtain  a  general  idea  of  the  validity  of  the  continuum 
equations  in  the  overall  flowfield  around  a  physical  object,  L  could  be  set  to  some 
characteristic  physical  dimension  of  the  object  (e.g.  the  chord  of  a  wing  on  an  aircraft). 
On  the  other  hand,  to  obtain  an  understanding  of  the  validity  of  the  continuum  equations  in 
approximating  the  fluid  physics  of  the  flow  field,  L  should  be  set  to  a  characteristic  gradient 
length  of  a  fluid  property,  e.g.  Figure  1.2  illustrates  generally  accepted  bounds  on  the 
regions  of  validity  for  the  continuum  equation  sets.  Note  that  the  Boltzmann  equation, 
the  governing  equation  of  kinetic  theory,  exhibits  no  Knudsen  number  dependence  for 
sufficiently  dilute  gases. 
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Figure  1.2:  Knudsen  Number  Validity  Ranges  for  Various  Equation  Sets  [65] 
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The  branch  of  gas  dynamics  which  is  formally  concerned  with  accounting  for  changes 
in  the  gas  which  occur  due  to  interactions  at  the  molecular  level  is  kinetic  theory.  Although 
the  basic  ideas  of  kinetic  theory  originated  in  the  mid-1800’s  with  Ludwig  Boltzmann  [23], 
the  practical  solution  and  application  of  these  methods  remain  areas  ripe  for  mathematical 
research.  The  complexity  involved  with  performing  a  kinetic  analysis  of  a  gas  has  presented 
a  challenging  problem  from  the  day  these  methods  were  conceived. 

Kinetic  theory  attempts  to  track,  at  least  statistically,  the  energy  state,  momenta  and 
position  of  every  particle  in  the  gas  as  a  function  of  time,  and  accounts  for  the  variation 
of  these  properties  due  to  collisions.  This  information  can  then  be  integrated  over  the 
collection  of  particles  to  obtain  the  macroscopic  properties  of  the  gas.  No  inherent 
assumptions  regarding  equilibrium  are  required  and  only  the  particle  collision  dynamics 
require  modeling.  The  assumptions  which  go  into  such  models  typically  do  not  exhibit 
restricted  validity  in  regions  of  non-equilibrium. 

Like  continuum  fluid  mechanics,  various  governing  equations  have  been  derived  to 
describe  the  behavior  of  molecules  in  a  gas  with  various  assumptions  regarding  interactions 
and  collision  dynamics.  The  most  popular  and  well  studied  of  these  is  the  Boltzmann 
equation  (alternatives  include  for  example  the  Kac  master  equation  [45]). 

Unfortunately,  the  great  flexibility  afforded  by  kinetic  theory  is  not  without  cost.  The 
Boltzmann  equation  is  a  nonlinear  integro-differential  equation  for  a  probability  density 
function  which  statistically  describes  the  energy  state  of  the  particles  as  a  function  of 
time.  For  a  monatomic  gas,  this  equation  must  be  solved  in  a  space  of  dimension  no  less 
than  seven.  For  polyatomic  molecules,  which  may  possess  several  components  of  angular 
momentum  and  vibrational  degrees  of  freedom,  the  dimension  of  the  space  grows  even 
higher. 

The  present  work  considers  a  gas  which  is  monatomic  and  further  assumes  that  the 
internal  energy  capacity  of  such  molecules  is  negligible.  In  other  words,  particles  will 
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be  assumed  to  possess  only  translational  energy.  This  restriction  allows  study  of  the 
phenomena  of  interest  without  adding  unnecessary  complexity.  Under  this  assumption, 
recall  that  the  translational  kinetic  energy  of  a  particle,  e,  is  related  to  its  velocity  via 

e  =  ^m\\v\\2  (1.3) 

Here,  ||  •  ||  represents  the  standard  Euclidean  norm.  In  light  of  (1.3),  one  may  equivalently 
examine  thermodynamic  properties  in  terms  of  the  distribution  of  molecular  velocities 
rather  than  energy. 

Consider  a  gas  occupying  some  physical  region  VQ  M3  during  some  interval  of  time. 
The  velocity  density  function,  f  :  V  x  M3  x  M+  — >  M+,  is  a  probability  density  function 
that  describes  the  distribution  of  velocity  over  the  collection  of  particles.  The  function  is 
defined  such  that 

/  (r,v,t)  drdv  (1.4) 

represents  the  probability  of  finding  a  particle  in  the  differential  element  of  physical  space, 
dr  centered  at  r,  which  possesses  velocity  in  the  differential  element  of  velocity  space,  dv 
centered  at  v,  at  time  t.  Here,  and  in  all  integrals  to  follow,  the  use  of  the  notation  dr  is  used 
to  mean  dr{dr2dr3  and  likewise  the  term  dv  is  used  for  dv\dv2dv3. 

It  is  at  this  point  that  one  first  may  begin  to  grasp  the  large  dimensionality  associated 
with  examining  a  gas  at  the  microscopic  level  as  /  is  seen  to  be  defined  over  a  seven¬ 
dimensional  space.  In  comparison,  the  variables  involved  in  the  partial  differential 
equations  modeling  the  gas  at  the  macroscopic  level  are  defined  at  most  over  a  four¬ 
dimensional  space.  The  six  dimensions  of  physical  and  velocity  space  are  commonly 
referred  to  as  the  “phase  space”.  It  should  also  be  noted  that  many  variations  of  /  appear 
in  the  literature  with  some  authors  referring  to  Nf  or  nf  as  the  distribution  function,  where 
N  is  the  number  of  molecules  in  the  gas  and  n  is  the  number  density  of  molecules  in  the 
gas.  For  clarity  in  this  work,  /  will  always  be  taken  as  the  probability  density  function  for 
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molecular  velocity,  or  what  some  authors  have  referred  to  as  the  normalized  distribution 
function. 

For  any  function,  Q ,  of  molecular  velocity,  v,  we  define  the  expectation  value  Q  at 
position  r  and  time  t  as 

Q(r,t)=  f  Q(v) f  (r,v,t)dv  (1.5) 

Such  functions  are  of  great  importance,  as  they  allow  for  the  computation  of  the 
macroscopic  thermodynamic  properties  of  the  gas.  This  is  illustrated  in  Table  1.1  [41]. 


Table  1.1:  Expectation  Values  Related  to 
Macroscopic  Properties. 


Q  (r,  v,  t) 

Q 

Description 

N 

n 

Number  Density 

mN 

P 

Mass  Density 

Nv 

n 

u 

Fluid  Velocity 

-**ln  [/] 

S 

Entropy 

3^'-  “I2 

T 

Temperature 

fllv-Ml2 

P 

Pressure 

One  can  therefore  see  the  importance  of  the  density  function  in  determining  the 
macroscopic  thermodynamic  properties  of  interest  in  the  fluid.  It  is  with  this  in  mind 
that  one  seeks  a  relation  to  describe  the  evolution  of  the  density  function  throughout 
phase  space.  One  such  description  is  provided  by  the  Boltzmann  equation  which  accounts 
for  changes  to  /  due  to  three  influences:  particle  advection,  acceleration  of  particles  by 
external  forces,  and  intermolecular  collisions.  The  equation  may  be  modified  to  include  the 
distribution  of  energy  over  various  internal  energy  modes,  but  for  simplicity  we  consider 
only  the  basic  case  of  a  simple,  monatomic  gas.  In  this  case,  the  Boltzmann  equation  is 
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given  by 


o 

■7~f  ( r ,  v,  t)  +  v  •  V,-/  (r,  v,t)  +  F  ■  Vv/  (r,  v,t)  =  J  [/]  (r,  v,  f) 


(1.6) 


Here,  F  is  any  externally  applied  force,  and  V,.  and  V,  represent  the  gradient  with 
respect  to  position  and  velocity  variables,  respectively.  That  is,  Vr  =  and 

Vv  =  ■  The  collision  integral  /  is  defined  as 

J[f](r,v,t)  =  f  f  [f(r,v'(v,w,n),t)f(r,w'(v,w,n),t)~ 

Jr3  Js+ 

f  (r,  v,  t)  f  ( r ,  vv,  0]  B  ( 6 ,  ||v  -  w||)  dddedw 

where  S+  denotes  the  positive  half  of  the  unit  sphere  in  M3,  n  is  the  unit  vector  defining  the 
collision  orientation,  6  and  e  are  the  azimuth  and  elevation  angles  of  n,  respectively,  and  B 
is  the  collision  cross  section.  The  post-collision  velocities  V ,  vv'  are  given  by 

v' (v,w,n)  =  v  -  [n  ■  (v  -  w)]n  (1.7) 

W  (v,  w,  n)  =  w  +  [n  •  (v  -  w)]  n  (1.8) 

Given  the  complexity  inherent  to  (1.6),  it  is  not  difficult  to  understand  why  closed  form 
solutions  of  the  Boltzmann  equation  for  flows  of  practical  interest  are  difficult  to  obtain. 
The  Boltzmann  equation  will  be  discussed  in  greater  detail  in  Chapter  2. 

1.2  Computational  Methods  for  the  Boltzmann  equation 

Computational  methods  for  the  Boltzmann  equation  continue  to  be  developed; 
however,  due  to  the  inherent  complexity  of  the  equation  itself,  such  methods  are 
correspondingly  complicated.  This  section  summarizes  the  strengths  and  weaknesses  of  a 
few  of  the  methods  currently  employed  in  this  field.  These  methods  can  largely  be  broken 
into  two  categories:  deterministic  methods  and  stochastic  methods. 

1.2.1  Deterministic  Methods. 

Deterministic  methods  for  the  Boltzmann  equation  are  typically  plagued  by  two  issues 
over  their  stochastic  counterparts,  namely,  increased  complexity  and  degraded  physical 
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fidelity.  On  the  other  hand,  they  avoid  the  introduction  of  variance  due  to  stochastic 
processes. 

1.2.1. 1  Direct  Boltzmann  CFD. 

A  class  of  methods  known  as  Direct  Boltzmann  solvers  or  Direct  Boltzmann 
Computational  Fluid  Dynamics  (CFD)  borrows  from  the  methods  of  continuum  CFD. 
These  methods  attempt  to  track  the  evolution  of  the  velocity  distribution  function  by 
performing  finite  difference  or  finite  volume  computations  over  a  grid  in  phase  space. 

A  number  of  issues  arise  with  such  an  approach.  Developing  a  numerical  grid  in  the 
six-dimensional  phase  space  is  much  more  complex  than  simply  gridding  physical  space 
as  in  continuum  CFD.  When  discretized  over  a  velocity  grid,  the  distribution  function  is 
constrained  to  have  compact  support.  In  general,  this  is  not  physically  accurate  and  it  is 
non-trivial  to  establish  appropriate  bounds  on  the  support  region  a  priori.  Furthermore,  grid 
refinement  in  velocity  space  is  non-trivial,  and  as  of  yet  an  area  of  open  research,  leading 
to  the  requirement  of  including  a  large  number  of  grid  points  in  order  to  obtain  an  accurate 
representation. 

Even  with  these  drawbacks,  these  methods  have  a  number  of  desirable  features  that  are 
not  necessarily  present  in  other  computational  schemes.  First,  they  are  directly  traceable 
to  the  Boltzmann  equation.  Second,  the  numerical  methods  employed,  although  tailored 
to  the  Boltzmann  equation,  are  well  understood  in  general  from  a  mathematical  viewpoint, 
and  are  directly  amenable  to  stability  and  error  analysis. 

Work  on  such  methods  began  in  the  late-1960’s  and  has  continued  to  the  present  [7, 
48].  Due  to  the  associated  complexities  and  computational  demands  required  for  accuracy, 
these  methods  have  not  yet  gained  popular  acceptance  as  practical  engineering  tools  but 
have  been  successfully  employed  to  solve  basic,  fundamental  flows. 
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1.2. 1.2  Molecular  Dynamics. 

Molecular  dynamics  methods  have  been  employed  since  the  late  1950’s  [6,  16]. 
The  idea  behind  such  methods  is  to  simulate  the  evolution  of  a  very  large  collection 
of  particles  throughout  the  flowfield.  Particle  interactions  and  collisions  are  computed 
deterministically. 

The  complexity  associated  with  such  methods  is  fairly  substantial.  First,  as  these 
methods  seek  to  deterministically  describe  the  evolution  of  the  particles,  each  particle’s 
trajectory  must  be  advanced  by  computing  all  of  the  possible  field  interactions  with  all  other 
particles.  Secondly,  as  collision  interactions  are  computed  deterministically,  selection  of  an 
appropriate  collision  partner  cannot  be  decoupled  from  particle  advection/field  interactions. 
Finally,  scatter  due  to  the  finite  sample  size  decreases  only  as  one  over  the  square  root  of 
the  sample  size  [16].  Therefore,  in  order  to  achieve  accuracy,  a  large  number  of  particles 
must  be  simulated. 

As  such,  molecular  dynamics  approaches  are  more  frequently  employed  when  the 
number  of  simulated  particles  can  be  nearly  on  the  order  of  the  number  of  actual  particles 
in  the  flow.  In  these  cases,  the  molecular  dynamics  method  yields  accurate  results  and  has 
been  applied  to  the  study  of  elementary  physical  problems. 

1.2. 1.3  Discretization  Methods. 

There  is  a  class  of  deterministic  methods  known  as  discretization  methods  which  seek 
to  “discretize”  the  behavior  of  the  gas  molecules.  One  such  method  is  known  as  the  Lattice 
Gas  Automata  Method.  In  this  method,  particle  advection  through  phase  space  is  viewed 
as  motion  from  one  node  to  an  adjacent  node  in  physical  and  velocity  space.  A  similar 
method  known  as  the  Discrete  Velocity  Method  discretizes  only  velocity  space.  That  is,  a 
finite  set  of  possible  particle  velocities  is  specified  [17]. 

Such  methods,  while  interesting  in  their  own  right,  cannot  fully  replicate  the  physics 
represented  in  the  Boltzmann  equation.  Namely,  it  is  impossible  to  totally  conserve  mass, 
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momentum,  and  energy  over  a  finite  velocity  set  [17].  This  leads  to  a  requirement  for 
inclusion  of  a  large  number  of  potential  velocities  and  the  computations  associated  with 
determining  the  most  appropriate  post  collision  velocities  over  such  a  set  is  quite  complex. 
Although  consistency  is  an  issue  for  such  methods  it  has  been  shown  that  when  a  large 
enough  set  of  possible  velocities  is  taken,  these  methods  can  be  accurate  for  elementary 
problems.  Their  associated  complexity  has,  however,  rendered  them  impractical  for 
engineering  applications  at  this  time  [38]. 

1.2.2  Stochastic  Methods. 

Stochastic  methods  for  the  Boltzmann  equation  typically  employ  a  simulation  algo¬ 
rithm  of  relatively  lower  complexity  than  deterministic  methods  to  model  intermolecular 
interactions.  For  this  reason,  initial  attempts  at  developing  such  methods  were  aimed  not  at 
preserving  mathematical  consistency  with  the  Boltzmann  equation  as  much  as  accurately 
simulating  the  physics  of  particle  interactions.  For  this  reason,  convergence  results  on  such 
methods  are  somewhat  limited.  Additionally,  the  incorporation  of  various  stochastic  pro¬ 
cesses  in  such  simulations  leads  to  the  introduction  of  variance  in  the  solution  that  is  not 
present  in  their  deterministic  counterparts. 

1.2.2. 1  The  Direct  Simulation  Monte  Carlo  Method. 

The  Direct  Simulation  Monte  Carlo  (DSMC)  method  was  originally  developed  in 
the  mid-1960’s  by  Bird  [15].  The  method  relies  on  the  Monte  Carlo  method  originally 
developed  by  von  Neumann  and  Ulam  (although  published  in  1950  by  Forsythe  and 
Liebler  [40]).  The  method  is  based  on  a  stochastic  simulation  of  a  fraction  of  the  actual 
number  of  particles  in  the  gas.  Each  simulated  particle  is  taken  to  represent  N/Np  actual 
particles. 

Unlike  the  molecular  dynamics  method,  particle  interactions  are  computed  probabilis¬ 
tically.  The  principal  approximation  is  termed  the  uncoupling  principle,  which  allows 
intermolecular  collisions  to  be  decoupled  from  particle  advection.  Collision  interactions 
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appropriate  for  the  time  increment  are  computed  probabilistically  after  which  particles  are 
advected  along  their  velocity  vectors.  Each  simulated  collision  represents  N/Np  actual  col¬ 
lisions.  Use  of  a  grid  in  physical  space  ensures  that  the  selected  collision  interactions  are 
between  near  neighbors  [17].  Specifically,  one  can  view  this  as  a  two-stage  splitting  scheme 
with  collisionless  advection  integrating 

df 

i+vVr/  =  0  (1.9) 

ot 

and  the  collision  simulation  integrating 

%=J[f](r,v,t)  (1.10) 

Initially,  the  DSMC  method  was  met  with  some  consternation.  Although  remaining 
true  to  the  principles  of  kinetic  theory,  the  method  itself  was  not  formally  derived  from  the 
Boltzmann  equation.  Without  a  proof  of  consistency  or  convergence,  this  raised  questions 
as  to  the  validity  of  the  method  itself.  By  the  early  1980’s,  however,  Nanbu  [58]  and  others 
had  proposed  methods  derived  from  the  Boltzmann  equation  itself,  and  by  1992  Bird’s 
method  had  been  proven  to  be  consistent  with  the  Boltzmann  equation  as  well  [75]. 

The  DSMC  method  is  not  without  its  drawbacks.  First,  a  significant  number  of 
particles  must  be  simulated  to  achieve  results  representative  of  reality.  This  raises 
storage  issues  as  the  position  and  velocity  vectors  of  each  simulated  particle  must 
be  stored  throughout  the  simulation.  Selection  of  representative  collisions  over  these 
potentially  large  data  arrays  introduces  a  significant  burden  on  the  simulation  process. 
Secondly,  the  stochastic  nature  of  the  simulation  introduces  a  significant  amount  of 
variation  in  the  results  and  in  practice  data  must  be  averaged  over  an  ensemble  to 
reduce  the  variance  in  the  solution.  Even  with  these  drawbacks,  DSMC  is  the  standard 
computational  method  employed  when  increased  accuracy  over  the  continuum  equation 
sets  is  required  and  has  gained  acceptance  in  the  field  of  hypersonic  aerodynamics  [24,  44, 
73],  flows  involving  micro  electromechanical  systems  (MEMS)  [37],  and  semiconductor 
manufacturing  processes  such  as  film  deposition  processes  [71]. 
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With  only  a  single  velocity  vector  per  simulated  particle,  DSMC  approximates  the 
velocity  density  function  utilizing  a  point  measure  approximation  as  illustrated  by 

1  Np 

f(y)  =  —  V  <s(v  -  v/)  (l.ii) 

1SP  i= 1 

where  v,-  is  the  velocity  vector  of  the  ith  particle. 

Equation  1.11  was  the  starting  point  for  Nanbu’s  derivation  of  an  algorithm  consistent 
with  the  Boltzmann  equation.  Within  each  cell  of  the  grid,  the  solution  is  assumed 
to  be  space  homogeneous.  Beginning  with  the  weak  form  of  the  space  homogeneous 
Boltzmann  equation,  Nanbu  was  able  to  develop  appropriate  stochastic  processes  that 
governed  selection  of  collision  pairs  and  collision  outcomes  within  a  given  cell.  Such 
methods  will  be  discussed  in  greater  detail  in  Chapter  3. 

1.2.2.2  Stochastic  Particle  Methods. 

Recently,  Wagner  sought  to  generalize  the  DSMC  method  into  a  class  of  methods  he 
termed  Stochastic  Particle  Methods  [61,  62,  76].  His  approach  was  motivated  by  a  desire 
to  allow  simulated  particles  to  have  a  varying  level  of  influence,  possibly  depending  on 
their  location  in  the  flowfield.  This  technique  would  potentially  improve  the  ability  of 
DSMC  to  handle  flowfields  which  contain  both  low  and  high  density  regions  of  simulated 
particles.  Specifically,  as  flow  expands  around  an  obstacle,  it  is  common  that  the  density 
of  simulated  particles  in  the  wake  region  fall  to  levels  that  exorbitantly  increase  the  level 
of  scatter  present  in  the  solution  [62].  This  is  predominantly  due  to  the  fact  that  although 
very  few  simulated  collisions  are  being  calculated,  they  will  have  tremendous  impact  on  the 
approximation  to  the  density  function,  since  each  simulated  collision  may  in  fact  represent 
millions  of  actual  collisions.  Secondly,  the  very  assumption  that  several  million  particles 
share  the  same  velocity  vector  in  a  region  where  very  few  simulated  particles  exist  leads 
to  a  gross  underrepresentation  of  the  density  function  itself.  While  the  stochastic  particle 
method  developed  by  Wagner  attempts  to  address  the  first  of  these  issues,  it  does  little  to 
combat  the  second,  as  the  method  continues  to  employ  point  measures.  This  is  a  goal  of 
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the  work  presented  herein  that  will  be  discussed  in  much  greater  detail  throughout  this 
document. 

Specifically,  the  stochastic  particle  method  approximates  the  density  function  as 
illustrated  by 

Np(t) 

f(v,t)=^gi(t)6(v-Vj)  (1.12) 

(=1 

where  gj  is  a  weighting  term  on  the  influence  of  the  z'-th  particle.  Also,  note  that  Np  is 
now  allowed  to  vary  with  time.  The  reason  for  this  is  due  to  another  feature  of  the  method, 
namely  a  collision  softening  effect  that  requires  Np  in  general  to  grow  as  more  collisions 
are  computed.  Specifically,  each  collision  generates  two  new  particles  which  are  assigned 
a  weight  and  post-collision  velocity,  while  the  two  original  particles  persist  with  their 
velocities  unaltered,  but  with  their  weighting  reduced.  The  post-collision  weightings  of  the 
four  particles  must  be  chosen  in  such  a  way  as  to  ensure  conservation  of  mass,  momentum, 
and  energy.  The  key  criticism  of  the  method  is  that  because  of  this,  the  number  of  simulated 
particles  grows  with  time  which  is  highly  undesirable  from  an  implementation  standpoint. 
To  combat  this,  a  method  was  proposed  by  which  particles  would  be  discarded  after  their 
weight  fell  below  a  certain  threshold  [61].  Nevertheless,  perhaps  because  of  the  complexity 
associated  with  the  creation  and  destruction  of  simulated  particles,  the  method  has  not  seen 
widespread  application. 

1.2.2. 3  Low  Variance  Deviational  Simulation  Monte  Carlo. 

Recently,  Baker,  Homolle,  and  Hadjiconstatinou  have  developed  a  simulation 
method  which  employs  Monte  Carlo  techniques  only  to  simulate  the  departure  from 
thermodynamic  equilibrium  [10,  43].  This  method  is  termed  Low  Variance  Deviational 
Simulation  Monte  Carlo  (LV-DSMC).  The  method  was  specifically  developed  to  handle 
low  Mach  number  flows  (e.g.  microscale  applications)  in  such  a  way  as  to  attempt  to 
ensure  the  variance  is  independent  of  the  magnitude  of  the  departure  from  thermodynamic 
equilibrium.  With  traditional  DSMC,  a  limiting  factor  is  that  for  near-equilibrium  solutions, 
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large  numbers  of  simulated  particles  must  be  employed  to  reduce  the  variance  to  the  point 
that  the  deviation  is  observable. 

This  technique  is  based  on  defining  the  deviation,  fd ,  from  an  arbitrary  Maxwell- 
Boltzmann  (equilibrium)  density  function  as 


f  =  f~fMB 


(1.13) 


where  / 'mb  is  a  Maxwell-Boltzmann  density  of  the  form  given  by 


,  ,  ,  hMb 

fMB  (v)  =  7^— v  exp 


/2 nkTuB\ 


m  (v  -  Umb) 


2  kT 


MB 


(1.14) 


where  nMB,  T mb ,  and  uMb  are  the  number  density,  temperature  and  mean  velocity  associated 
with  the  Maxwell-Boltzmann  density  function.  As  shown  in  any  text  on  kinetic  theory, 
the  Maxwell-Boltzmann  density  is  the  solution  of  the  Boltzmann  equation  in  translational 
equilibrium.  Substituting  (1.13)  into  the  (1.6),  a  new  form  of  the  collision  integral  is 
obtained.  Baker,  Homolle,  and  Hadjiconstantinou  were  then  able  to  develop  a  particle 
scheme  with  particle  advection  and  collision  substeps  to  simulate  the  transformed  equation. 

The  process  is  more  complicated  than  traditional  DSMC,  as  particles  do  not 
necessarily  advect  according  to  (1.9),  and  their  advection  behavior  is  dependent  upon  how 
fMB  is  chosen.  As  for  the  collision  substep,  it  is  viewed  as  having  a  two  step  effect:  first, 
a  net  change  to  fMB,  and  second,  a  change  to  fd  with  the  goal  being  to  lump  as  much  of 
the  change  as  possible  into  fMB  and  then  regenerate  deviational  particles  with  velocities 
corresponding  to  the  new  fd.  Therefore,  like  the  stochastic  particle  method,  LVDSMC  also 
relies  on  particle  creation  and  destruction  at  each  time  step. 

The  summary  of  methods  provided  here  is  in  no  way  totally  comprehensive.  In  fact, 
many  methods  currently  being  developed  blur  the  lines  of  distinction  outlined  here  [16] 
(e.g.  Monte  Carlo  based  Discrete  Velocity  Methods  [56]).  Each  method  presents  unique 
challenges  in  its  implementation  and  representation  of  the  Boltzmann  solution.  In  this 
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work,  we  focus  on  developing  a  stochastic  particle  based  method  which  borrows  from  the 
ideas  of  DSMC,  but  removes  the  point  measure  representation  of  the  density  function. 

1.3  Overview  of  Current  Approach  -  Distributional  Monte  Carlo 

This  work  outlines  the  development  of  a  new  approach  for  the  solution  of  the 
Boltzmann  equation  which  has  been  termed  “Distributional  Monte  Carlo”.  This  method 
falls  under  the  category  of  stochastic  particle  methods,  but  unlike  attempts  by  prior  authors, 
the  Distributional  Monte  Carlo  method  removes  the  point  measure  approximation  of  the 
density  function  by  allowing  each  simulated  particle  to  possess  not  just  a  single  velocity 
vector,  but  rather,  a  complete  velocity  distribution  function.  Binary  collisions  computed 
between  simulated  particles  have  the  effect  of  altering  the  particles’  density  function. 
Specifically,  a  binary  collision  between  two  simulated  particles  is  viewed  as  a  space 
homogeneous  relaxation  of  the  distribution  function  of  the  potentially  millions  of  actual 
particles  represented  by  the  two  simulated  particles;  rather  than  assuming  that  millions  of 
collisions  each  have  the  same  outcome  (as  in  traditional  DSMC). 

The  reasons  for  such  modifications  are  many.  First,  the  assumption  that  the  millions 
of  actual  particles  represented  by  a  single  simulated  particle  all  possess  the  same  velocity 
vector,  is  nonphysical.  Kinetic  theory  tells  us  that  intermolecular  collisions  will  drive 
the  velocity  density  function  of  the  collection  towards  the  Maxwell-Boltzmann  density 
function.  Second,  in  addition  to  being  nonphysical,  this  assumption  has  the  effect  of 
over  emphasizing  collision  effects  and  introducing  additional  variance  in  the  solution 
particularly  in  areas  where  the  density  of  simulated  particles  is  low  as  observed  by 
Wagner  [76].  Third,  the  point  measure  approximation  limits  convergence  to  occur  only 
in  the  weak  sense  to  the  Boltzmann  solution. 

To  illustrate  the  restrictive  nature  of  the  point  measure  representation,  consider  an 
arbitrary  collision  pair  of  simulated  particles.  Figure  1.3(a)  is  a  2-dimensional  conceptual 
illustration  of  the  two  simulated  particles  as  well  as  the  collection  of  actual  particles 
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represented  by  each  simulated  particle.  Note  that  every  particle  represented  by  a  given 
simulated  particle  possesses  the  same  velocity  vector.  Figure  1.3(b)  illustrates  the  velocity 
density  function  of  the  entire  collection  of  actual  particles  prior  to  collision.  The  function 
is  singular  exactly  at  the  velocities  of  the  two  simulated  particles  and  identically  zero 
everywhere  else. 


(a)  A  Collision  Pair  Prior  to  Collision 
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(b)  Velocity  Density  Function  of  Actual  Particles 
Figure  1.3:  DSMC  Pre-Collision  Modeling 
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Figure  1.4(a)  illustrates  the  same  collision  pair  post  collision.  Note  that  the  collision 
process  changes  the  velocity  vectors  of  all  the  particles  and  all  of  the  actual  particles 
represented  by  a  simulated  particle  are  assumed  to  possess  the  same  velocity  vector  post 
collision.  Figure  1 .4(b)  illustrates  the  velocity  density  function  of  the  collection  of  actual 
particles  post-collision.  The  function  is  still  singular  at  the  new  post-collision  velocities 
and  identically  zero  elsewhere.  The  locus  of  possible  collision  outcomes  is  given  by  the 
overlaid  circle.  That  is  to  say,  any  two  velocities  which  lie  on  the  circle  directly  across  from 
one  another  represent  a  possible  collision  outcome.  As  there  are  infinitely  many  possible 
outcomes,  the  selection  of  a  single  outcome  is  performed  stochastically. 


The  Distributional  Monte  Carlo  Methods  developed  by  the  author  represent  the  density 
function  not  as  in  ( 1 . 1 1 )  or  ( 1 . 1 2),  but  rather  as 

1  Np 

/(v)  =  — Vy;(v)  (i.i5) 

iyP  ;=  l 

where  f  is  the  density  function  of  the  ith  simulated  particle.  In  developing  such  a  scheme 
for  the  space  homogeneous  case,  the  most  important  choices  one  must  consider  are  how 
to  approximate  f  and  how  to  compute  collision  interactions  between  particle  pairs.  The 
current  work  presents  the  results  for  two  specific  schemes. 

The  first  scheme  is  an  approach  which  employs  a  fixed  functional  form  of  f,  and  is 
termed  DMC-KDE,  where  KDE  stands  for  Kernel  Density  Estimator.  In  this  case,  each 
particle’s  density  function  assumes  a  prescribed  form,  namely 


where,  K  e  Lr  (M3)  satisfies 


K  (v)  dv  =  1 

3 

vK  (v)  dv  =  0 


(1.16) 


(1.17) 

(1.18) 
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(a)  A  Collision  Pair  Prior  to  Collision 


Post-Collision  Velocity  Density  Function 


(b)  Velocity  Density  Function  of  Actual  Particles 

Figure  1.4:  DSMC  Post-Collision  Modeling 


and  h  6  M+.  Under  these  conditions,  (1.15)  assumes  the  form  of  a  Kernel  Density 
Estimator,  where  K  is  the  kernel  function  and  h  is  the  bandwidth.  Unlike  a  point 
measure  approximation,  it  will  be  shown  that  an  approach  based  on  this  technique  achieves 
convergence  for  L°°  (m3)  solutions  of  the  space  homogeneous  Boltzmann  equation,  as  well 
as  pointwise  convergence  for  bounded  solutions,  whereas  DSMC  converges  only  in  law  (or 
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weakly).  A  sufficient  criteria  for  convergence  of  the  method  will  be  shown  that  h  must  be 
chosen  as  a  function  of  Np  in  such  a  way  that  h  —>  0  as  Np  — »  oo. 

An  intuitive  choice  for  K  is  a  Maxwellian  density,  in  which  case  the  actual  particles 
represented  by  a  single  particle  are  no  longer  assumed  to  possess  the  exact  same  velocity, 
but  rather  have  velocities  distributed  according  to  a  Maxwellian  distribution.  The  physical 
interpretation  of  this  approach  is  that  the  actual  particles  represented  by  a  simulated 
particle  are  in  thermodynamic  equilibrium  with  one  another.  This  situation  is  illustrated 
in  Figure  1.5(a).  Note  that  the  actual  particles  are  no  longer  constrained  to  possess  the 
same  velocities.  The  velocity  density  function  prior  to  collision  of  the  actual  particles 
represented  by  the  simulated  collision  pair  is  illustrated  in  Figure  1.5(b).  Note  that  unlike 
DSMC,  the  distribution  function  is  no  longer  singular,  but  rather  bimodal,  being  a  sum  of 
Maxwellians. 

As  will  be  discussed  later,  the  DMC-KDE  method  represents  a  kernel  density 
estimator  applied  to  the  DSMC  technique.  While  simulated  particles  are  allowed  to 
posses  non-singular  distributions,  collision  simulation  is  still  performed  using  similar  rules 
to  DSMC  to  assign  the  mean  of  the  post-collision  Maxwellian  distributions.  The  post¬ 
collision  situation  is  outlined  in  Figure  1.6(a).  In  this  case  the  actual  particles  are  not 
assumed  to  possess  the  exact  same  velocity,  but  note  that  the  mean  velocity  of  the  collection 
is  the  same  as  the  post-collision  velocity  computed  in  DSMC  (Figure  1.4(a)).  The  post¬ 
collision  velocity  density  function  is  given  in  Figure  1.6(b).  In  this  case,  valid  solutions  for 
the  center  points  of  the  Maxwellian  density  functions  will  lie  on  opposite  sides  of  the  locus 
of  the  overlaid  circle. 

As  will  be  discussed,  the  replacement  of  the  point  measure  representation  with  the 
non-singular  DMC-KDE  form  allowed  the  author  to  prove  strong  convergence  of  the 
method  for  L°°  (M3)  and  bounded  solutions  of  the  space  homogeneous  Boltzmann  equation. 
This  was  the  first  time  that  a  stochastic  particle  simulation  was  shown  to  converge  in  this 
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(a)  A  Collision  Pair  Prior  to  Collision 

Pre-Collision  Velocity  Density  Function 

0.04 

0.035 

0.03 

0.025 

0.02 

0.015 

0.01 

0.005 

0 

-3-2-10123 

vx 

(b)  Velocity  Density  Function  of  Actual  Particles 

Figure  1.5:  DMC-KDE  Pre-Collision  Modeling 


sense,  as  opposed  to  simply  in  law.  Unfortunately,  it  will  be  observed  that  the  method  does 
not  achieve  a  significant  variance  reduction  over  traditional  DSMC.  Even  so,  the  method  is 
valuable  from  a  practical  perspective  in  that  it  allows  for  direct  evaluation  of  the  velocity 
density  function. 

Although  valuable  in  its  own  right,  DMC-KDE  does  not  fully  embody  the  Distribu¬ 
tional  Monte  Carlo  concept  envisioned  by  this  research.  In  the  general  Distributional  Monte 
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(a)  A  Collision  Pair  Post  Collision 


Post-Collision  Velocity  Density  Function 


(b)  Velocity  Density  Function  of  Actual  Particles 

Figure  1.6:  DMC-KDE  Post-Collision  Modeling 

Carlo  approach,  particles  may  possess  arbitrary  velocity  density  functions,  and  interactions 
are  computed  as  a  space  homogeneous  relaxation  over  the  current  time  step  of  the  com¬ 
bined  density  function  of  the  two  simulated  particles  involved  in  a  collision.  Figures  1.7(a) 
and  1 .7(b)  illustrate  the  pre-collision  conditions  for  two  simulated  particles  each  possessing 
Maxwellian  distribution  functions  prior  to  collision. 
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(a)  A  Collision  Pair  Pre  Collision 
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(b)  Velocity  Density  Function  of  Actual  Particles 

Figure  1 .7 :  DMC  Pre-Collision  Modeling 


Instead  of  basing  the  collision  modeling  on  a  binary  collision  between  the  simulated 
particles,  the  method  computes  the  post-collision  velocity  density  function  of  both  particles 
by  computing  an  approximate  space  homogeneous  relaxation  of  the  initial  density  function 
over  the  current  time  step.  As  the  time  step  is  increased  the  post  collision  solution  tends 
toward  a  Maxwellian.  This  is  illustrated  in  Figures  1.8(a)  and  1.8(b). 
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(a)  A  Collision  Pair  Post  Collision 
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(b)  Velocity  Density  Function  of  Actual  Particles 

Figure  1.8:  DMC  Post-Collision  Modeling 


In  this  work,  the  author  derives  and  proves  that  so  long  as  the  binary  collision 
relaxation  calculation  is  consistent  with  the  Boltzmann  equation,  the  Distributional  Monte 
Carlo  method  converges  in  law  to  the  Boltzmann  solution.  Additionally,  the  method 
exhibits  strong  convergence  for  L°°  (R3)  solutions.  This  marks  the  first  development  of 
a  non-point  measure  based  stochastic  particle  method  for  the  Boltzmann  equation,  as  well 


25 


as  the  first  proof  of  such  a  method’s  convergence.  This  is  the  main  contribution  of  this 
work. 

Finally,  in  order  to  obtain  numerical  results,  the  method  was  implemented  using  a 
BGK  [11]  approximation  for  collision  simulation.  This  approach  is  termed  DMC-BGK. 
Although  not  completely  consistent  with  the  Boltzmann  equation,  the  BGK  approximation 
is  commonly  applied  in  rarefied  gas  dynamics.  It  should  be  noted  that  this  does  not  detract 
from  the  theoretical  proof  previously  discussed,  as  many  suitable  approaches  exist  which 
are  consistent  with  the  space-homogeneous  Boltzmann  equation  (e.g.  moment  methods, 
spectral  methods,  etc.);  the  BGK  approximation  was  chosen  simply  because  its  ease  of 
coding  allowed  for  rapid  generation  of  numerical  solutions.  The  approach  is  applied 
to  the  well-studied  Bobylev-Krook-Wu  [22,  50]  problem,  where  it  is  observed  to  have 
significantly  reduced  variance  over  traditional  DSMC. 

The  remainder  of  this  work  is  outlined  as  follows:  Chapters  2  and  3  provide  the 
necessary  background  and  results  on  the  space-homogeneous  Boltzmann  equation  and  the 
Direct  Simulation  Monte  Carlo  Method.  Chapter  4  presents  the  DMC-KDE  method,  proof 
of  convergence  of  the  method,  and  numerical  results  on  the  Bobylev-Krook-Wu  problem. 
Chapter  5  presents  the  general  Distributional  Monte  Carlo  approach  in  detail  as  well  as 
proof  of  its  convergence  in  law  to  the  Boltzmann  solution.  Also  in  Chapter  5,  the  DMC- 
BGK  approach  is  detailed  and  applied  numerically  to  the  Bobylev-Krook-Wu  problem. 
Chapter  6  summarizes  the  conclusions  and  findings  of  this  work  and  outlines  potential 
areas  of  fertile  research  for  future  investigators. 
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II.  The  Boltzmann  equation 


The  Boltzmann  equation  (1.6)  is  the  governing  equation  of  kinetic  theory.  As 
introduced  in  Chapter  1,  the  fundamental  assumption  of  kinetic  theory  is  that  all 
macroscopic  properties  of  a  gas  can  be  deduced  from  a  knowledge  of  the  interactions 
and  internal  structure  of  its  constituent  molecules  [41].  Remarkably,  the  fundamentals  of 
kinetic  theory  were  first  formally  postulated  by  Ludwig  Boltzmann  [23]  at  a  time  when 
the  atomic  makeup  of  matter  was  not  an  accepted  concept.  Modern  physics  realizes 
Boltzmann’s  vision  as  the  application  of  statistical  mechanics  to  a  gas.  Although  the 
development  of  these  concepts  can  be  traced  back  nearly  150  years,  it  has  been  only 
over  the  last  50  years  that  significant  contributions  have  been  made  to  understanding 
some  of  the  mathematical  properties  (e.g.  existence  and  uniqueness  of  solutions)  of  the 
Boltzmann  equation  [1,  2,  26,  46,  54].  Numerical  methods  for  the  Boltzmann  equation 
were  successfully  developed  and  applied  to  the  study  of  basic  physical  flows  over  the  last 
30  years  [14,  18,  57],  but  have  only  recently  made  significant  inroads  to  applications  of 
practical  interest  in  both  hypersonic/high  altitude  flight  [24,  44]  and  microscale  flows  [37]. 

To  understand  and  appreciate  the  development  of  the  methods  described  in  this  work, 
it  is  necessary  to  begin  with  a  brief  general  discussion  of  the  Boltzmann  equation  (1.6).  This 
chapter  briefly  outlines  the  assumptions  behind  and  derivation  of  the  Boltzmann  equation, 
before  simplifying  to  the  space  homogeneous  Boltzmann  equation  which  is  the  basis  for 
the  remainder  of  this  work.  The  chapter  concludes  with  the  presentation  of  a  well  known 
solution  [22,  50]  of  the  space  homogeneous  Boltzmann  equation  which  has  been  employed 
in  numerous  studies  [52,  64]  to  examine  the  performance  of  various  numerical  methods. 

The  Boltzmann  equation  is  a  probabilistic  representation  of  the  evolution  of  the  energy 
state  of  the  molecules  comprising  the  gas  under  consideration.  In  general,  a  molecule  may 
store  its  energy  in  a  number  of  modes.  It  may  possess  energy  in  translational  and  rotational 
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motion.  It  may  also  store  energy  in  a  vibrational  mode  along  one  or  more  chemical  bonds. 
Additionally,  the  molecule  may  become  electronically  excited.  To  simplify  the  analysis 
of  the  present  work,  it  will  henceforth  be  assumed  that  the  gas  under  consideration  is 
monatomic  and  further,  that  the  internal  energy  capacity  of  such  molecules  is  negligible.  In 
other  words,  particles  will  be  assumed  to  possess  only  translational  energy.  This  restriction 
allows  study  of  all  of  the  phenomena  of  interest  without  adding  unnecessary  complexity 
and  does  not  prohibit  later  generalization  to  molecules  with  internal  structure. 

As  the  Boltzmann  equation  is  a  probabilistic  model  of  the  gas,  some  basic  concepts 
from  probability  theory  must  be  introduced  before  continuing  with  a  brief  derivation  of  the 
Boltzmann  equation.  The  following  definitions  motivate  the  concept  of  a  probability  space. 

Definition  2.1.  Given  a  set  Q.  and  a  cr-algebra  S  of  subsets  of  X,  a  countably  additive 
function,  p,  from  S  into  [0,  oo]  is  called  a  measure  and  (Q,S,p)  is  called  a  measure  space. 


Definition  2.2.  A  measurable  space  (Q,  S)  is  a  set  Q  with  a  cr-algebra  S  of  subsets  ofQ. 

Definition  2.3.  Given  a  measurable  space  (Q,  S),  a  probability  measure,  P,  is  a  measure 
on  S  with  P  ( Q)  =  1.  (Q,S,P)  is  called  a  probability  space.  Elements  of  S  are  called 
events. 

Definition  2.4.  A  law  is  a  Borel  probability  measure;  that  is  a  probability  measure  defined 
on  the  Borel  cr-algebra. 

The  Boltzmann  equation  is  a  relation  describing  the  evolution  of  a  law  on  the 
measurable  space  (D,D),  where  D  =  V  x  K3,  and  D  is  the  Borel  cr-algebra  defined  on  D. 
Here,  V  c  M3  is  the  physical  volume  occupied  by  the  gas,  and  the  remaining  dimensions 
of  D  represent  three  components  of  velocity.  In  rarefied  gas  dynamics,  D  is  commonly 
referred  to  as  the  phase  space.  The  following  theorem  is  necessary  to  motivate  the  concept 
of  probability  density. 
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Theorem  2.1.  (Radon-Nikodym)  [35]  On  the  measurable  space  (X,  S)  let  p  be  a  cr-finite 
measure.  Let  v  be  a  finite  measure,  absolutely  continuous  with  respect  to  p.  Then  for  some 


he£l  C X,S,p ), 


v(E) 


£ 


hdp 


(2.1) 


for  all  E  G  S.  Any  two  such  h  are  equal  p-almost  everywhere,  h  is  called  the  Radon- 


Nikodym  derivative  or  density  ofv  with  respect  to  p,  and  is  written  h  =  E. 


Definition  2.5.  A  law  P  on  Rk  is  said  to  have  a  density  f  if  P  is  absolutely  continuous 
with  respect  to  Lebesgue  measure  Ak  and  has  Radon-Nikodym  derivative  =  f.  In  other 
words,  P  (A)  =  jAf(x)  dAk  ( x)for  all  Borel  sets  A. 


For  purposes  of  this  work,  the  term  probability  density  will  be  used  to  describe  the 
Radon-Nikodym  derivative  as  given  in  Definition  2.5  so  as  to  avoid  confusion  with  the 
physical  properties  of  the  gas  of  mass  density  or  number  density.  As  mentioned  above,  the 
Boltzmann  equation  describes  the  evolution  of  a  probability  law  for  molecular  velocity.  To 
be  more  precise,  the  Boltzmann  equation  is  an  integro-differential  equation  over  the  space 
D  for  the  evolution  of  the  probability  density  of  the  law  for  molecular  velocity.  Physically 
speaking,  the  probability  of  a  molecule  existing  in  a  given  subset  of  phase  space  A  c  D  is 
given  by  P  (A)  =  fAf(x)  dx,  where  dx  denotes  the  Lebesgue  measure  in  the  phase  space. 
This  is  equivalent  to  the  more  common  differential  description  that  /  (r,  v)  drdv  represents 
the  probability  of  molecule  existing  in  the  differential  element  dr  centered  at  r  in  physical 
space,  while  possessing  velocity  in  the  differential  element  dv  centered  at  v  in  velocity 
space.  Having  established  the  required  mathematical  definitions,  a  brief  derivation  of  the 
Boltzmann  equation  based  on  physical  reasoning  is  next  presented. 


2.1  Overview  of  Derivation  of  the  Boltzmann  Equation 

A  brief  exposition  of  the  origin  of  the  Boltzmann  equation  and  the  development 
of  relations  that  will  later  prove  useful  is  discussed  herein.  Detailed  derivations  and 
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discussions  of  the  Boltzmann  equation  can  be  found  in  numerous  sources  (see,  for  example 
Reference  [23,  27-30,  32,  41,  74]). 

There  are  a  few  important  physical  assumptions  that  are  inherent  to  the  Boltzmann 
equation.  First,  one  must  assume  complete  collisions.  That  is,  the  “dt”  time  interval 
involved  in  the  equation  is  much  larger  than  the  collision  interaction  time  [41].  This 
is  somewhat  more  precise  and  restrictive  than  the  more  common  interpretation  that  only 
binary  collisions  occur  in  the  gas.  The  assumption  implies,  in  fact,  that  the  collisions  occur 
so  rapidly  (instantaneously)  that  the  probability  of  a  tertiary  collision  is  zero  [41]. 

The  next  assumption  is  that  /  is  slowly  varying  at  small  scales,  i.e.  on  the  order  of 
molecular  diameters [41].  This  assumption  is  based  on  the  concept  of  a  thermodynamic 
limit,  namely,  that  at  the  physical  scales  of  interest  the  gas  must  be  capable  of  being 
described  in  a  statistically  meaningful  way.  Unless  the  gas  is  quite  dense  this  assumption 
has  no  implication  on  how  /  varies  on  the  order  of  the  mean  free  path.  This  assumption 
limits  the  Boltzmann  equation  to  dilute  gases. 

The  final  assumption  inherent  to  the  Boltzmann  equation  is  the  existence  of  molecular 
chaos  [41].  This  means  that  no  correlation  can  exist  between  the  molecular  velocities  of 
any  two  particles  which  are  outside  of  one  another’s  field  of  influence.  These  assumptions 
can  be  shown  to  be  consequential  from  a  single  assumption  termed  the  Boltzmann- Grad 
limit  [29,  30]. 

The  Boltzmann  equation  accounts  for  changes  to  the  density  function  throughout  the 
phase  space  by  three  mechanisms:  particle  advection  through  physical  space,  particle 
acceleration  by  means  of  some  external  force,  and  intermolecular  collisions.  Let  drdv 
represent  a  differential  element  of  phase  space  centered  at  (r,  v).  By  accounting  for  these 
mechanisms,  one  can  determine  a  balance  equation  for  particle  accumulation  within  the 
differential  element  of  phase  space  over  the  differential  time  dt.  Heuristically,  we  seek  an 
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expression  of  the  form 


df 

dt 


O,  v,  t) 


■m. 


(r,  v,  t)  + 


df\ 
dt  J 


(r,  v,  t)  + 


Accel 


df\ 
dt l 


0,  v,  t ) 


(2.2) 


Co// 


Accel 


where,  (|Q  represents  the  rate  of  change  of  /  due  to  particle  advection, 
represents  the  rate  of  change  of  /  due  to  particle  acceleration,  and  ^  represents  the  rate 

of  change  of  /  due  to  intermolecular  collisions.  Considering  only  elementary  mechanics,  it 
can  be  shown  [29,  74]  that  the  changes  due  to  advection,  and  particle  acceleration  are  given 
by 


{ ( r V’  0  =  _V  •  Vr./  (r>  V,  0  (2.3) 

l  ^  )  Com' 

(r,  v,  t)  =  -F  (r,  t)  ■  Vv/  (r,  v,  t)  (2.4) 

I  Accel 


where  F  (r,  /)  is  the  total  external  force  acting  at  r  at  time  t. 

The  term  resulting  from  collisions  requires  some  additional  considerations.  At  a  given 
time  and  location  in  phase  space,  we  must  account  for  all  possible  collisions  that  result 
in  particles  with  pre-collision  velocities  outside  the  differential  element  dv  centered  at  v 
attaining  a  post-collision  velocity  within  dv.  Additionally,  we  must  account  for  all  possible 
collisions  that  result  in  particles  with  pre-collision  velocities  within  the  differential  element 
dv  centered  at  v  attaining  post-collision  velocities  outside  dv.  Generically,  we  will  represent 
these  gain  and  loss  terms  by  the  symbols  G  and  L,  respectively,  and  write 


df\ 

dt) 


(r,  v,  t)  =  G  ( r ,  v,t)-L  (r,  v,  t) 


(2.5) 


Coll 


In  order  to  analyze  G  and  L  properly,  one  must  first  consider  the  mechanics  of  a  binary 
collision  between  two  particles.  Conceptually,  the  easiest  model  to  employ  is  the  hard 
sphere,  or  billiard  ball  model  of  a  molecule.  Although  the  model  is  physically  imperfect,  it 
provides  a  credible  starting  point  for  this  analysis  and  the  modifications  required  for  other 
molecular  models  will  be  discussed  shortly. 
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To  begin,  consider  two  particles  undergoing  a  collision  with  velocities  v  and  w 
respectively.  Denote  the  relative  velocity  of  the  particles  by  V  =  v-w.  Denote  the  velocities 
of  the  two  molecules  post  collision  by  v'  and  w' ,  and  the  post  collision  relative  velocity  by 
V'  =  v'  -  w'. 

During  the  collision,  both  momentum  and  kinetic  energy  must  be  conserved,  which 
implies  [30] 


v  +  w  =  V  +  W  (2.6) 

INI2  +  INI2  =  Hv'll2  +  llw'U2  (2.7) 

Let  n  be  a  unit  vector  directed  along  v  -V .  Then  the  post-collision  velocities  are  related  to 
the  pre-collision  velocities  via 


V  =  v  -  [n  •  V]  n 

(2.8) 

w'  =  w  +  [«•!/]  n 

(2.9) 

which  implies 


V'  =  V-  [2/7  •  V]  n 


(2.10) 


From  the  reference  frame  of  the  first  molecule,  this  means  that  the  second  molecule  is 
specularly  reflected  and  n  is  directed  along  the  line  of  centers.  This  reference  frame 
provides  a  convenient  method  for  computing  the  gain  and  loss  terms  sought  above. 
Specifically,  let  d  represent  the  diameter  of  a  single  molecule.  In  this  reference  frame 
let  particle  1  be  a  sphere  at  rest  and  endowed  with  diameter  2d,  and  view  the  remaining 
particles  as  point  masses  moving  with  velocity  V  =  w  -  v.  Next,  one  must  determine 
the  probability  that  a  particle  possessing  velocity  w  impacts  the  sphere  at  the  differential 
surface  element  dS  =  d2dn  during  the  time  interval  [/,  t  +  dt\.  If  a  particle  is  to  exhibit  such 
a  collision  it  must  lie  within  the  slant  cylinder  traced  out  by  the  differential  element  over  dt 
which  has  height  ||V]|  dt  and  volume  \V  ■  n\  dSdt.  The  probability  that  a  molecule  satisfies 
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these  requirements  is  given  by 


f{r,  w,  t) \V  •  n\  dSdwdt  (2.11) 

where  dw  represents  a  differential  element  of  velocity  space  centered  at  w. 

What  we  have  established  in  (2.1 1)  is  the  probability  of  a  collision  between  a  molecule 
with  velocity  in  the  differential  element  dw  impacting  a  single  molecule  with  velocity  in 
the  differential  element  dv,  but  have  said  nothing  of  the  probability  of  the  first  molecule 
possessing  a  velocity  in  the  differential  element  dv.  Therefore  the  probability  of  any 
molecule  with  velocity  in  the  differential  element  dw  impacting  any  molecule  with  velocity 
in  the  differential  element  dv  per  unit  time  is  given  by  [74] 

/(r,  v,  t)  f  {r,  w,t)\V  ■  n\  dSdwdv  (2.12) 

This  is  commonly  referred  to  as  a  collision  of  class  (v,  w)  [74].  The  L  term  we  seek  in  (2.5) 
is  due  to  the  loss  of  particles  with  velocity  in  the  differential  element  dv  due  to  collisions 
with  particles  of  any  other  velocity.  Therefore,  we  conclude 

L(r,  v,  t)  =  I  |  f  (r,  v,  t)f(r,  w,  t)  d2 \V  ■  n\  dndw  (2.13) 

Jr3  Js- 

where  S~  is  the  hemisphere  corresponding  to  V  ■  n  <  0,  i.e.  orientations  over  which 
particles  are  moving  toward  each  other  prior  to  the  collision  [29] .  A  similar  analysis  can 
be  performed  to  determine  an  expression  for  G.  To  do  this,  one  must  consider  the  concept 
of  inverse  collisions  or  replenishing  collisions  for  a  given  collision  class.  Given  a  collision 
between  molecules  with  velocities  Vi  and  V2,  an  inverse  collision  is  any  collision  which  has 
initial  velocities  equal  to  the  final  velocities  of  the  original  collision  and  the  same  direction 
of  the  line  of  centers  [74].  Using  this  concept  one  can  determine  the  G  term  we  seek  in 
(2.5),  which  represents  the  replenishment  of  particles  in  the  differential  element  dv  centered 
at  v  due  to  collisions,  is  given  by  [30,  74] 

G(r,v,t)=  I  I  f  (r,  v' ,t)  f  (r,w' ,t) d2 \V  ■  n\dndw  (2.14) 

Jr3  Js+ 
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Therefore,  we  conclude  that 

[%)  ( r,v,t)=  f  f  [/ (r,v',t)f  (r,w',t)  — 

\ dt  I cou  Jr3 

f  (r,  v,  t)  f  (r,  w,  0]  d 2  |(v  -  w )  •  n\  dndw 

where, 

V  -v  -  [n  ■  (v  -  w)]  n 
w'  -  w  +  [n  •  (v  -  w)]  n 


(2.15) 


(2.16) 

(2.17) 


The  right  hand  side  of  (2.15)  is  known  as  the  collision  integral  and  is  the  main  source  of 
complexity  in  the  Boltzmann  equation.  Thus,  the  Boltzmann  equation  for  a  monatomic 
hard  sphere  gas  is  given  by 

^  (r,  v,  t)  +  v  •  V,./  (r,  v,  t)  +  F  (r,  t)  ■  Vv/  (r,  v,  t)  = 

(2.18) 


dt 


f  f  [f(r,v',t)f(r,w',t)-f(r,v,t)f(r,w,t)\d2\(v-w)-n\dndw 

Jr3  Js+ 

As  previously  mentioned,  other  molecular  models  beyond  hard  sphere  exist;  some  of 
which  are  preferred  due  to  improved  physical  accuracy  and  others  because  of  mathematical 
simplicity.  It  can  be  shown  [28,  29]  that  in  such  cases  the  only  required  modification  to 
the  above  is  to  replace  the  term  d2  |(v  -  w)  ■  n\  with  a  function  of  ||v  -  w||,  and  6,  the  angle 
between  n  and  y  -  w.  In  this  case,  the  Boltzmann  equation  becomes 


(r,  v,t)  +  v-  Vr/  (r,  v,t)  +  F  (r,  t )  •  V,,/  (r,  v,  t)  = 
dt 


(2.19) 


f  f  [f(r,v',t)f(r,w',t)-f(r,v,t)f(r,w,t)\B(0,\\v-w\\)dwdOde 

JR3  JS  + 

where  e  is  the  second  angle  defining  the  orientation  of  n  on  the  unit  sphere,  and  B 
is  called  the  collision  cross  section  or  collision  kernel.  In  the  case  of  hard  spheres, 
B  ( 6 ,  ||v  -  vv'j|)  =  cos  ( 6 )  sin  ( 6 )  ||v  -  w\\.  Another  common  model  is  called  the  inverse  power 
law,  in  which  the  attractive  force  between  two  particles  is  assumed  to  vary  as  the  n-th 
inverse  power  of  the  distance  between  them.  In  this  case,  B  is  given  by 


B  (6,  ||v  —  w||)  =  [3  (69  1 1 v  -  w||» 


(2.20) 
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where  the  function  (3  incorporates  the  collision  orientation.  The  case  of/?  =  5  gives  rise  to 
the  Maxwell  molecule,  in  which  case  the  dependence  on  ||v  -  w\\  vanishes  [53]. 

2.1.1  Relation  to  the  Continuum  Equations  of  Fluid  Mechanics. 

During  this  general  discussion  of  the  Boltzmann  equation,  it  is  worth  noting  its 
connection  to  the  continuum  equations  of  fluid  motion. 

Upon  taking  the  moment  of  the  Boltzmann  equation  with  respect  to  m,  mv,  and  \m\\v\\2, 
expressions  for  the  conservation  of  mass,  momentum,  and  energy  in  the  gas  will  result. 
These  expressions  will  vary  with  respect  to  the  form  of  /.  If  /  is  locally  Maxwellian,  the 
resulting  set  of  conservation  equations  reduces  to  the  well  known  Euler  equations  of  fluid 
mechanics. 


%  +  V-pu  =  0 
dt 

(2.21) 

+  (V  •  pu)  u  +  VE  =  0 

(2.22) 

+  V  ■  (e  +  P)u  =  0 
dt 

(2.23) 

where  e  is  the  specific  energy.  As  stated  previously,  however,  an  equilibrium  (locally 
Maxwellian)  description  is  not  sufficient  to  capture  all  of  the  phenomena  of  interest.  A 
method  which  attempts  to  address  this  situation  was  proposed  first  by  Hilbert  [42]  and 
later  developed  independently  by  Chapman  and  Enskog  [31,  36].  The  method  expands  the 
distribution  function  in  a  Hilbert  series  from  a  local  Maxwellian. 


/  =  /m(i+®(1)  +  0(2)  +  ---) 


(2.24) 


Truncating  the  expansion,  one  may  develop  expressions  for  the  O(0  terms.  For 
example,  including  the  first  correction  term  one  obtains  the  Navier-Stokes  equations. 

dp 


dt 


+  V  •  pu  =  0 


dpu  .  __  1  _ 

—  +  (V  •  pu)  u  +  Vp  =  — VTjj 
dt  Re 

de  1 

_  +  v.(e  +  p)„  =  -_[v.«] 


(2.25) 

(2.26) 
(2.27) 
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The  inclusion  of  the  first  correction  term  allows  for  the  physical  phenomena  of 
viscosity  and  thermal  conductivity  in  the  gas.  This  results  in  the  appearance  of  the  shear 
stress  tensor,  r,7,  the  heat  flux  vector,  q,  as  well  as  the  Reynold’s  Number,  Re,  and  Prandtl 
Number,  Pr.  Including  the  second  order  term  yields  the  Burnett  equations,  the  third  yields 
the  Super-Burnett  equations  and  so  forth.  It  is  important  to  note  however  that  while  the 
connection  is  seen  to  exist  to  the  continuum  equation  sets,  there  is  no  guarantee  that  such 
a  truncated  approximation  will  be  a  good  one.  In  fact,  open  debate  continues  over  whether 
the  Burnett  and  higher  order  sets  add  any  accuracy  over  the  Navier-Stokes  equations, 
and  it  has  been  shown  that  there  are  situations  in  which  they  violate  the  second  law  of 
thermodynamics  [34,  39]. 


2.2  The  Space  Homogeneous  Boltzmann  equation 

Clearly,  much  of  the  complexity  associated  with  the  Boltzmann  equation  is  inherent  to 
the  right-hand  side  which  is  known  as  the  collision  integral.  For  this  reason,  it  is  common 
to  study  the  space  homogeneous  version  of  the  equation.  It  is  also  common  to  neglect  the 
influence  of  external  forces,  in  which  case  the  equation  reduces  to 

d~f(v,t)=  f  f  [f(v',t)f(w',t)-f(v,t)f(w,t)\B(6,\\v-w\\)dwd6de  (2.28) 

Ot  Jr3  Js  + 

The  equation  is  solved  subject  to  the  initial  condition  f(v,  0)  =  /0  (v),  where  /0  e  L1  (M3) 
is  a  non-negative  function  which  satisfies  the  normalization  condition 


dv  =  1 


(2.29) 


Additionally,  in  order  to  maintain  finite  energy,  the  boundary  condition, 


lim  /(v,  t)  =  0  for  all  t  >0  (2.30) 

IldH00 

is  imposed.  In  the  space  homogeneous  case,  any  changes  to  /  are  solely  due  to 
intermolecular  collisions.  Existence  and  uniqueness  of  solutions  of  (2.28)  was  established 
by  Akeryd  [1,2]. 
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Another  common  notation  is  to  denote  the  collision  integral  by  Q  (/,  /)  where  Q  is  the 
functional  defined  by  [30] 


Q(f,g)=  f  f  [f(y',t)g(w',t)  +  f(w',t)g(V,t) 

Jr3  Js+ 

-f  (v,  t)  g  (w,  t)-  f  (w,  t)  g  (v,  t )]  B  (6, 1 1 v  -  w||)  dwdOde 
In  this  case  the  space  homogeneous  equation  can  be  written  as 

ot 


(2.31) 


(2.32) 


A  fundamental  solution  is  the  case  when  Q  (/,  /)  =  0.  In  this  case,  the  density  function 
is  constant  and  depleting  collisions  are  exactly  balanced  by  replenishing  collisions  for  all 
collision  classes.  This  state  is  known  as  translational  thermodynamic  equilibrium.  It  can  be 
shown  [74]  that  the  solution  for  which  this  criteria  holds  is  the  Maxwell-Boltzmann  density 


function 


/(v)  =  o-i  exp  [~a2  ||v  -  m||2] 


(2.33) 


where  <*1,2  are  constants  that  depend  upon  the  physical  characteristics  of  the  gas,  and  u 
represents  the  bulk  fluid  velocity.  One  frequently  considers  the  case  where  u  -  0  which  is 
known  as  a  nondrifting  Maxwellian  [29] . 

Let  (p  be  a  real  valued  function  over  M3  such  that  JT,  (p  (v)  Q  (/,  /)  dv  exists.  Cercignani 


shows  [29] 


<P  (v)  Q  (/,  /)  dv  =  - 


[f  W)  f  {W)  -  f  (y)  f  (w)] 


(2.34) 


[cp  (v)  +  (p  (w)  -  cp  (v')  -  <p  (w')]  B  (6,  ||v  -  w||)  dwdvdn 


as  well  as  the  more  compact  form 


(p  (v)  Q  (/,  /)  dv  = 


/(v)/(w)' 


(2.35) 


[cp  (v')  +  (p  (’ w ')  -  (p  (v)  -  cp  (w)]  B  ( 6 ,  ||v  -  w||)  dwdvdn 
The  right-hand  side  of  (2.35),  can  be  made  to  be  zero  independent  of  /  if 


0  (v)  +  (p  (w)  =  (p  (v')  +  (p  (wr) 


(2.36) 
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almost  everywhere  in  velocity  space  [29].  Since  the  left-hand  side  of  (2.35)  represents 
the  net  rate  of  change  of  the  average  value  of  (p  due  to  collisions,  any  (p  for  which 
J^3  (p  (v)  Q  (/,  /)  dv  =  0  for  any  /  is  termed  a  collision  invariant. 

It  can  be  shown  [29]  that  there  exists  a  five-dimensional  subspace  of  continuous 
functions  which  are  collision  invariants  that  is  spanned  by  the  functions,  {<A/}f=1 ,  given  by 


<Ai  (v)  =  1 

(2.37) 

^2  (V)  =  Vx 

(2.38) 

<MV)  =  Vy 

(2.39) 

<Mv)  = 

(2.40) 

<Mv)  =  llvll2 

(2.41) 

This  is  not  a  surprising  result,  as  the  appearance  of  </q  is  a  statement  of  the  conservation 
of  mass,  whereas  ^4)  result  from  the  conservation  of  momentum  in  each  of  the 

cardinal  directions  and  ip5  is  a  statement  of  the  conservation  of  energy. 

Next,  we  note  that  a  special  case  exists  when  (p  =  log  /.  If  Jr3  log  (/)  Q  (/,  /)  dv  exists, 
it  can  be  shown  [74]  that 


Q(f,f)dv<  0 


(2.42) 


This  is  known  as  the  Boltzmann  inequality  and  is  required  to  prove  the  Boltzmann  H- 
Theorem,  which  applies  not  just  to  the  space  homogeneous  equation  but  to  the  general 
Boltzmann  equation  as  well. 


Theorem  2.2.  Boltzmann ’s  H-Theorem.  Let  'IT  (t)  =  fR3f  (v,  t)  log  (/  (v,  t ))  dv  Then 

d<H 

—  (t)  <  0  (2.43) 

dt 

Further,  ^  =  0  if  and  only  if  f  is  a  Maxwell- Boltzmann  density  function. 


The  function  Fi  is  commonly  referred  to  as  Boltzmann’s  H-function.  This  theorem 
is  the  analogue  of  the  second  law  of  thermodynamics  at  the  microscopic  level,  in  fact, 


38 


it  can  be  shown  [47,  74]  that  the  entropy  of  the  gas  is  related  to  the  H-function  by 
S  =  This  is  important  because  it  demonstrates  that  the  Boltzmann  equation 

possesses  the  basic  thermodynamic  feature  of  irreversibility  [30].  It  also  indicates  that 
the  Maxwellian  distribution  is  the  prevailing  distribution  in  thermodynamic  equilibrium,  a 
state  characterized  by  zero  entropy  generation. 

2.2.1  Results  on  Existence  and  Uniqueness  of  Solutions. 

Existence  and  uniqueness  of  solutions  for  the  space  homogeneous  Boltzmann  equation 
(2.28)  is  a  well  studied  area.  Carleman  [25]  proved  existence  and  uniqueness  of  solutions 
for  continuous  initial  conditions  which  satisfied  (1  +  ||v|0/o  (v)  6  L°°  (M3)  for  some  k  >  6. 
Morgenstem  [54,  55]  proved  the  same  for  Maxwellian  molecules  with  a  truncated  collision 
kernel  and  initial  conditions  /0  e  Ll  Povzner  [59]  proved  existence  and  uniqueness 
for  continuous  collision  kernels  that  obey  B{6,  V)  <  C(1  +  V),  where  C  e  1'  and  initial 
conditions  /0  with  (1  +  ||v||4)/0(v)  6  L1  (M3).  The  results  presented  herein  are  due  to 
Akeryd  [1,  2],  and  represent  a  more  general  result  for  bounded  collision  kernels.  The 
following  theorem  combines  several  of  Akeryd’s  results. 


Theorem  2.3.  Suppose  that 

0  <  B(6,\\v-w\\)  <Kk  (2.44) 

for  all  v,w  6  M3,  6  e  [0,2^].  Here,  Kk e  M+.  If  fo  e  L1  (M3),  with  /0  >  0  and 
Jr3  Jo  (v)  dv  =  1,  then  there  exists  a  unique,  nonnegative  solution  f{v,t )  6  Ll  (R3)  of  the 
space  homogeneous  Boltzmann  equation  (2.28)  with  /(v,  0)  =  /0  (v). 

Furthermore,  if  (l  +  ||v||2)/o(v)  e  L1  (if3),  then  f(v,  t)(l  +  ||v||2)  6  Ll  for  all 
t  >  0  and 


r 

Jr3 


/  (v,  t)  dv 


r 

Jr 3 


vf  (v,  t )  dv  = 


'  f  (v,  t)  dv 


f 


fo  (v)  dv 

3 

vfo  (v)  dv 


Hv||2/o(v)  dv 


(2.45) 

(2.46) 

(2.47) 
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for  all  t  >  0.  Also,  iff0  log  /0  e  L1  (M3),  then  f  (v,  t )  log  /  (v,  0  e  L1  (U3  )/or  c//Z  t  >  0,  rmJ 

'K(r)  =  f  /  (v,  r)  log  (/  (v,  t))  dv  (2.48) 

J  R3 


is  a  nonincreasing  decreasing  function  in  t. 


2.2.2  The  Bobylev,  Krook,  and  Wu  Solution. 

Although  many  of  the  fundamental  mathematical  properties  such  as  existence  and 
uniqueness  of  solutions  of  the  Boltzmann  equation  have  been  solved,  and  several 
approximate  approaches  exist,  a  characterizing  feature  of  the  equation  has  been  the 
nearly  complete  absence  of  exact  analytical  results  [22].  One  well  known  exception  is 
the  Bobylev-Krook-Wu  solution  of  the  space  homogeneous  equation,  first  published  by 
Bobylev  [19-21]  and  later  independently  by  Krook  and  Wu  [49,  50].  The  problem  has  thus 
served  as  a  useful  test  case  for  various  schemes  (e.g.  [52,  58])  and  will  be  utilized  in  the 
same  capacity  by  the  work  presented  herein.  For  this  reason,  a  brief  background  on  the 
solution  is  presented  here. 

The  solution  is  specific  to  Maxwellian  Molecules.  Without  loss  of  generality,  Bobylev 
considered  the  following  initial  conditions 


r 

Jr3 


v/o  (v)  dv  -  0 


\\v\\2  fo(v)dv  =  3 


(2.49) 

(2.50) 


which  by  (2.41)  will  hold  for  all  t  >  0.  These  moments  correspond  to  the  net  momentum 
and  energy  of  the  gas,  and  thus  by  utilizing  (2.33)  the  corresponding  Maxwell-Boltzmann 
density  function  (the  solution  as  t  — >  oo)  can  immediately  be  determined  to  be 

/mb(v)  =  (2/rr3/2exp|--^J  (2.51) 


Bobylev’s  approach  was  to  perform  a  Fourier  transform  on  the  space  homogeneous 
Boltzmann  equation.  Under  this  representation,  the  velocity  density  function  is  replaced 
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by  the  Fourier  variable  ®  defined  as, 


®  (k,t)  =  f  f  (v,  t)  exp  (—ik  ■  v)  dv  (2.52) 

Jm3 

where  ke  M3  is  the  variable  in  the  Fourier  domain.  The  transformed  space  homogeneous 
equation  becomes 


<9® 

dt 


L 

L 


Q  (/,  /)  exp  (-ik  •  v )  dv 


P 


k  ■  n 

w 


® 


k  +  ||k||  n\  (k  —  ||k||  ?7 


—  ®  (0)  ®  (k) 


dn 


(2.53) 


Note  the  drastic  simplification  that  results  from  taking  the  Fourier  transform,  namely  the 
five  dimensional  integral  in  (2.28)  is  reduced  to  a  two  dimensional  integral.  The  strategy, 
therefore,  is  to  determine  the  solution  of  (2.53)  and  invert  the  transformation  to  obtain  the 
density  function  via  the  inversion  formula 


f(v,t)  = - -  f  ®  (k,  t)  exp  (ik  ■  v)  dk  (2.54) 

(2;r)3  Jr3 

(2.53)  must  be  solved  subject  to  the  transformed  initial  and  normalization  conditions 


®o(*) 


/o  (v)  exp  (-ik  ■  v)  dv 


(2.55) 


Also,  the  transformed  Maxwellian,  ®mb,  can  be  determined  by  direct  substitution  of  (2.51) 
into  (2.52)  and  is  given  by 


®mb  (k)  =  exp 


(2.56) 


Bobylev  then  sought  to  determine  a  similarity  solution  of  (2.53)  of  the  form 


®  (k,  t)  =  ®0  (ke  '“)  exp 


\\k\\ 


(l  -e-2"') 


(2.57) 


where  jue  M3.  This  is  further  simplified  by  seeking  an  isotropic  solution,  the  simplest  of 
which  can  be  written  as 


®  (x,  t)  =  (l  -  Qxe  exp  [-x  (l  -  0c  /l')J 


(2.58) 
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where  x  =  ^y-  and  the  parameter  0  is  any  real  number  in  the  interval  [ 0,  |J.  Note  the 
dependence  on  x  versus  k  clearly  exhibits  the  isotropic  property  of  the  solution.  The 
variable  Ae  R+  is  dependent  upon  the  collision  cross  section  and  is  given  by 


T  =  ^J  yS(6»)(l-cos  {Of)  d6  (2.59) 

Applying  the  inverse  transform  to  (2.58)  results  in  the  final  expression  for  the  density 
function 


/  (v, t )  = 


1 

'  v2  ' 

[  1  -Tit) 

flMI2 

_l\] 

(l/rr(0r3/2  P 

l  2rWJ 

[  T  W 

l  2r 

v\ 

(2.60) 


where  r(t)  =  1  -  Qe~A>.  This  was  the  first  known  nontrivial  closed  form  solution  of  the 
space  homogeneous  Boltzmann  equation  [22].  Bobylev  defines  the  following  normalized 
even  moments  of  the  distribution  function 


zn(t)  =  -  *  f  ||v||2”  f  (v,t)dv  n  =  1,2,...  (2.61) 

(2/7  +  1)!!  JR3 

where, 

n 

(2n+  1)!!  =  Y\(2i+  1)  (2.62) 

(=i 

In  the  process  of  obtaining  his  solution,  Bobylev  was  able  to  exploit  the  isotropic  nature 
of  the  solution  and  determine  expressions  for  these  moments  fairly  simply  in  the  Fourier 
domain.  The  resulting  normalized  moments,  z,n  are  given  by 

Zn  (0  =  (l  -  Qe~M)n~'  [l  +  in  -  1)  ®e-At]  (2.63) 


The  simplicity  of  this  closed  form  solution  lends  it  to  great  use  in  evaluating  numerical 
schemes  for  the  Boltzmann  equation  (e.g.  [52,  56,  58,  63,  64]),  and  will  be  used  towards 
that  end  in  the  present  work  as  well.  The  present  work  will  utilize  the  solution  for  A  =  1/6 
and  0  =  2/5.  Visualizing  the  velocity  density  function  can  be  difficult  given  its  inherent 
dimensionality,  however,  as  the  present  solution  is  spherically  symmetric  one  can  readily 
derive  the  density  function  for  molecular  speed  by  moving  to  spherical  polar  coordinates 
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and  integrating  out  the  directional  dependence  from  (2.60) 


t)  = 


An 


( 2n ) 


-3/2 


exp 


2t  it) 


1  + 


1  -  r 


(t)( 


T(t)  \2r  (t)  2 


(2.64) 


The  speed  density  function  is  shown  in  Figure  2.1  below.  In  many  situations  it  is  not 


Figure  2.1:  Molecular  Speed  Density  Function  of  Bobylev-Krook-Wu  Solution 

practical  to  directly  compare  the  density  function  (e.g.  DSMC,  which  converges  in  law  but 
does  not  directly  produce  a  convergent  density  function),  and  therefore  it  is  common  to 
employ  a  comparison  of  the  even  moments  of  the  Bobylev-Krook-Wu  solution.  The  first 
four  moments  are  illustrated  in  Figure  2.2. 

These  results  have  been  employed  in  the  current  work  as  a  numerical  test  case  for  the 
various  algorithms  under  development  and  represent  a  standard  benchmark  problem  for 
Boltzmann  simulators  and  solvers. 
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Figure  2.2:  First  Four  Even  Moments  of  Bobylev-Krook-Wu  Solution 
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III.  The  Direct  Simulation  Monte  Carlo  Method 


3.1  Overview 

The  Direct  Simulation  Monte  Carlo  Method  was  originally  proposed  by  Bird  in 
1963  [14].  Bird  entitled  the  method  Direct  Simulation  to  distinguish  from  methods  which 
would  be  characterized  as  Direct  Solution  techniques.  In  other  words,  Bird’s  original 
development  of  DSMC  was  never  intended  to  be  a  numerical  solution  technique  for 
the  Boltzmann  equation.  Bird  argued  that  a  computational  simulation  of  the  physics  of 
rarefied  gas  interactions  without  explicitly  incorporating  the  Boltzmann  equation  was  more 
tractable  and  perhaps  more  valuable  than  applying  numerical  techniques  to  the  equation 
itself  [17], 

Others  were  not  as  content  to  accept  the  simulation  process  without  a  connection  to 
the  underlying  equation  and  in  1980  Nanbu  proposed  a  DSMC  scheme  directly  derived 
from  the  Boltzmann  equation  [58].  By  1987,  Babovsky  [8]  had  proven  that  Nanbu’s 
scheme  was  in  fact  convergent  in  law  to  the  space  homogeneous  Boltzmann  solution.  Two 
years  later,  Babovsky  and  Illner  proved  the  same  for  the  nonhomogeneous  equation  [9]. 
Interestingly,  as  Bird’s  method  itself  argues  from  the  same  physical  assumptions  (e.g.  dilute 
gas,  molecular  chaos,  etc.)  inherent  to  the  Boltzmann  equation,  Wagner  was  later  able  to 
prove  that  in  spite  of  the  procedure’s  ignorance  of  the  Boltzmann  equation  itself,  it  did  in 
fact  converge  to  the  solution  of  the  Boltzmann  equation  in  law  [75]. 

We  utilize  results  which  stem  from  Nanbu’s  technique  and  Babovsky’s  proof  of  its 
convergence,  both  of  which  are  presented  in  this  chapter.  The  approaches  could  be  applied 
equally  well  to  develop  techniques  that  are  more  like  Bird’s  approach,  but  the  direct 
traceability  to  the  Boltzmann  equation  of  Nanbu’s  technique  results  in  a  relatively  easier 
path  to  proving  convergence.  For  this  reason,  we  utilize  results  from  Nanbu’s  technique 
over  those  of  Bird’s. 


45 


3.2  Nanbu’s  Simulation  Technique 

In  this  section,  an  overview  of  the  original  derivation  [58]  of  Nanbu’s  method  is 
provided  while  in  Section  3.3  an  explanation  [8]  is  given  which  is  more  amenable  to 
analyzing  convergence  of  the  method.  Although  other  DSMC  schemes  exist  which  are 
in  more  widespread  usage,  Nanbu’s  method  is  presented  here  due  to  the  clarity  with  which 
it  can  be  derived  from  the  Boltzmann  equation. 

The  basic  premise  of  any  DSMC  scheme  is  the  uncoupling  principle,  namely,  that 
particle  convection  and  particle  collisions  are  treated  independently.  In  non-homogeneous 
problems,  physical  space  is  discretized  by  means  of  a  volume  grid,  and  within  a  cell  the 
solution  is  assumed  to  be  space  homogeneous.  For  this  reason,  the  original  derivation 
of  Nanbu’s  scheme  focused  on  developing  the  stochastic  relations  to  simulate  the  space 
homogeneous  Boltzmann  equation.  Simulated  particle  convection  is  a  more  trivial 
problem,  as  particles  simply  convect  along  their  velocity  vectors  for  the  current  time  step 
by  A r  =  vA t.  The  space-homogeneous  simulation  of  the  collision  operator  is  therefore 
the  central  problem  in  DSMC.  As  such,  proof  of  convergence  of  Nanbu’s  method  was 
originally  obtained  for  the  space  homogeneous  equation  [8]  and  subsequently  for  the  full 
Boltzmann  equation  [9].  With  this  in  mind,  the  current  work  focuses  solely  on  space 
homogeneous  results. 

Nanbu’s  derivation  seeks  to  develop  a  set  of  stochastic  rules  governing  collision  pair 
selection  and  collision  outcomes  for  a  set  of  N  particles  each  possessing  a  velocity  vector 
V/.  Nanbu  begins  by  writing  the  initial  density  function  for  molecular  velocity  as  a  point 
measure  approximation 

^w  =  y;I>(v-'9  <3-‘> 

where  the  superscript  indicates  the  kth  time  step  of  the  simulation.  Expanding  the 
temporal  derivative  in  the  space  homogeneous  Boltzmann  equation  using  a  forward  Euler 
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discretization  in  time  and  solving  for  fk+l  yields 


fk+'  (v)  =  fk(v)  +  0  (/*,/*)  Ar 


Upon  substituting  (3.1)  into  (3.2),  Nanbu  [58]  shows  that  the  collision  integral  is  given  by 


where, 


Np  i=  i  j=  i 


Sjj=  |  #(v'  -  -  v*)  ||v  -  w||5(||v  -  w||  ,6)dndw  (3.4) 

JM3  js  + 

rh£l  8  (v  -  vk)  d(w  -  vk£j  ||v  -  w\\  2?(||v  -  w\\ ,  6)  dndw  (3.5) 

which  represent  the  replenishing  and  depleting  collision  contributions,  respectively.  Note 
that  Tkj  has  a  somewhat  simpler  form  than  S k-  and  may  be  evaluated  directly, 


r.HK-v‘ll45(v-v0 


where  Bk.  is  sometimes  referred  to  as  the  total  collision  cross  section  and  is  given  by 


In  order  to  simplify  Sku,  Nanbu  resorts  to  approximating  the  delta  function  as 


8  (v)  =  lim - —  exp 


Upon  substituting  (3.8)  and  (2.9)  into  (3.4),  one  obtains 


S  f,  =  lim - 

1  {ne) 


77  L 


GUv,  w) 


exP  [  e  Hvll2  +  H^ll2  + 


Kir + ikii-  -  (v- + v5)  ■  (v + w)  |  iiv  _  M\  dw 


where  Gk.  is  defined  by 


Gkj  (v,  w)  =  J'  exp  (akj  •  77)  B 


(llv  -  w\\,6)dn 


(3.10) 
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and  a\j  -  ^  ||v  —  w||  (v*  —  v*J.  As  stated  previously,  6  is  the  angle  between  (v  -  w )  and 
77.  Next,  introduce  the  spherical-polar  coordinate  system  with  the  axis  directed  along  ak-. 
Denote  the  angle  between  ak-  and  n  by  oj  and  the  azimuth  angle  of  n  about  the  new  axis  by 
Then  dn  =  sin  {oj)  dcod £,  and  6  =  6(co ,£).  Substituting  into  (3.10) 


X2  7r  r*n 

J  exp  [||a-;||  cos  (cu)]  B  (||v  -w\\,6  {oj,  £))  sin  {oj)  dojd% 


Nanbu  [58]  next  asserts  that  since  ak 


(3.11) 

oo  as  e  — >  0,  the  range  of  small  oj  is  the  dominant 
contributor  in  (3.11).  This  allows  him  to  utilize  the  following  approximations 

exp  [HU  cos  (m)]  *  exp  ||4||  |l  -  yj  (3.12) 

e(a>,&*0(  0,&=x  (3.13) 

sin(m)~u7  (3.14) 


whereof  is  the  angle  between  (v  -  w)  and  a*..  Nanbu  then  substitutes  these  approximations 
into  (3.1 1),  along  with  one  additional  approximation,  obtained  by  changing  the  upper  limit 
of  integration  from  n  to  oo  on  doj  to  obtain 


•2  71  r*°° 


Gku  (v,  w) 


r ^  r 
Jo  Jo 


exp 


a  1 


OJ 


ojB(\\v  -  w\\  ,x)  du)d% 


2nB(\\v  -w\\,x) 


exp 

( 4 ) 

ak. 

Uij 

Substituting  (3.16)  into  (3.9),  and  using  the  fact  that  ||a| 


li  =  lim 


f 

J  R3 


exp  (-4 


,J  e-»0  (ne) 

B  (||v  -  w\\  ,£)dw 

Recalling  (3.8),  this  expression  can  be  rewritten  as 

8 


(v  +  w  -  v)  -  vf)~  +  (| | v  -  w\\  -  || vt  -  v y 1 1 ) 


,  one  obtains 
21 


Sk  - 

O  ij  - 


vk  -  vk 

l  J 


J'  5  (v  +  w  -  Vi  -  6  (||v  -  w\\ 


|v/  —  v y ||)  B  (| | v  —  w| |  ,x)  dw 


vk  -  vk 

l  J 


vij 


vk  - 
*  J 


IK --111 -4) 


(3.15) 

(3.16) 


(3.17) 


(3.18) 
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where 


v*  +  v) 

v*k  =  v  -  - - J- 

<J  0 


(3.19) 


and xkij  is  the  angle  between  {yk  -  vk^j  and  v*k.  Substituting  (3.9),  (3.6),  and  (3.3),  into  (3.2) 
yields  the  final  expression  for  fk+1 


(3.20) 


p  i=  i 


where 


A  NP 

A  t 


P;  =  ~ 


N, 


P  7=1 


N„ 


k  4A t  vA 

Gf  =  fttZ<s 


P  7=1  V 


*  A: 

- 


vf-v5|h5(||vf-v;||,4) 


k  k 

v  —  V 

l  J 

(3.21) 

(3.22) 


Rewriting  Pk  as 


where, 


np 

-  V 

1=1 


P*  =  —  llv*  -  vA1|  fiA 

U  yy  II  '  7 II  IJ 


(3.23) 


(3.24) 


Nanbu  bases  his  scheme  on  the  interpretation  that  Pk  represents  the  probability  that  the 
ith  simulated  particle  undergoes  a  collision  during  the  time  interval  [4,4+ 1],  while  Pk- 
represents  the  probability  that  particle  i  collides  with  particle  j  during  At.  Under  this 
interpretation,  the  j  =  i  terms  must  be  omitted  from  the  summation,  as  a  simulated  particle 
cannot  collide  with  itself. 

Assuming  that  particle  i  collides  with  particle  j,  Nanbu  then  shows  that  the  conditional 
probability  density  function  for  collision  angle,  g,  is  given  by 

B(\\vk-vk\\,X)  sin  Or) 


g(x) 


Bk 

ij 


(3.25) 


Note  that  the  absence  of  4  from  the  expression  implies  that  the  azimuthal  dependence  is  of 
uniform  probability. 
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Based  on  these  results,  Nanbu  developed  the  following  stochastic  scheme  to  evolve 
the  molecular  velocities  during  a  given  time  step: 

•  For  each  particle,  calculate  Pk.  Generate  a  random  number  r\  in  the  interval  (0, 1) .  If 
Pk  >  r1?  accept  the  particle  for  collision. 

pk 

•  Sample  a  collision  partner  j  from  the  conditional  probability  distribution  P*m  =  -pf- , 

1  i 

by  sampling  a  second  random  number  r2  uniformly  from  the  interval  (0, 1)  and 
identifying  the  j  which  satisfies  X/J,  P*m  <  r2  <  X)Ui  P*m 


•  Sample  the  direction  of  v*  based  on  (3.25),  and  compute  the  post  collision  velocity 
of  the  ith  particle,  According  to 


v!  =  ^(llvf~v5ll*  +  vf  +  v5) 


(3.26) 


where  the  vector  R  is  a  unit  vector  computed  by  sampling  the  densities  for  azimuthal 
and  elevation  angles,  namely 

I  =  2tzt3  (3.27) 

where  r3  e  (0, 1)  is  another  random  number,^  is  determined  by  generating  a  random 
number  r4  e  (0, 1)  until  g  (£)  >  r4  where  x  —  nr  a,  and  R  is  computed  as 

i(f)cos  ( x ) 


R 


sin 

sin 


(?)  sin  Cl) 

(?) 


COS 


(3.28) 


In  general  Pk  ±  p1^  but  Nanbu  points  out  that  in  the  special  case  of  Maxwell 
molecules,  the  collision  probabilities  are  all  equal,  namely 

,(%bY(Np-r 


Pi  =  rfl  I 


\  m 


Nn 


A  Ni 


(3.29) 


where  yS0  and  b  are  parameters  which  define  the  collision  cross  section.  In  this  case, 
collision  partners  can  be  chosen  randomly  simplifying  the  second  step  of  the  procedure 
significantly. 
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Note,  however,  that  the  procedure  developed  by  Nanbu  does  not  exactly  produce  the 
solution  at  time  step  k+  1  from  equation  3.20.  Specifically,  note  that  the  effect  of  collisions, 
namely  the  Q1-  term,  does  not  lend  itself  to  be  written  as  a  simple  delta  function  involving 
two  distinct  post  collision  velocities.  Rather  the  form  of  this  term  indicates  that  the  post 
collision  velocities  for  a  collision  pair  are  distributed  over  a  sphere.  DSMC  techniques  such 
as  Nanbu’s  method  simply  sample  discrete  post  collision  velocities  from  this  sphere.  In  so 
doing,  much  of  the  information  contained  in  equation  3.20  is  destroyed.  If  one  considers  the 
fact  that  each  simulated  particle  represents  a  large  number  of  actual  particles,  the  solution 
indicates  that  the  collision  interaction  will  distribute  the  velocities  of  the  actual  particles 
over  the  sphere  post  collision.  Further,  the  appearance  of  the  first  term  in  equation  3.20 
indicates  that  a  portion  of  the  molecules  are  unaffected  by  the  collision  interaction.  As 
DSMC  allows  for  only  a  single  velocity  per  simulated  particle,  the  sampling  cannot  account 
for  all  of  these  effects. 

3.3  Proof  of  Convergence  of  Nanbu’s  Method 

The  previous  section  provided  an  overview  of  Nanbu’s  method  which  followed 
Nanbu’s  original  derivation  fairly  closely  [58].  This  section  will  follow  Babovsky’s  proof, 
which  showed  Nanbu’s  method  converged  in  law  to  the  solution  of  the  space  homogeneous 
Boltzmann  equation  (2.28)  [8].  Much  of  the  details  of  Babovsky’s  work  will  be  included 
here,  as  the  work  serves  as  a  good  template  for  completing  such  proofs,  and  the  following 
section  lays  out  many  of  the  details  omitted  from  Babovsky’s  original  paper.  Before 
continuing  we  define  this  form  of  convergence  concretely  as  follows. 

Definition  3.1.  Let  (S ,  T )  be  a  topological  space,  with  Borel  cr-algebra  B.  Define 
Cb(S )  :=  Cb(S,T )  to  be  the  set  of  all  bounded,  continuous,  real-valued  functions  on  S .  A 
sequence  of  laws  { Pn }  is  said  to  converge  to  a  law  P  if  fs  fdPn  — >  js  fdP  as  n  — »  oo  for 
every  f  in  Cb  (S)  [35 ]. 
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To  begin,  it  is  necessary  to  develop  a  few  additional  relations  and  consider  Nanbu’s 
approach  in  a  somewhat  more  rigorous  setting.  Recalling  the  expressions  for  the  post¬ 
collision  velocities  (2.9),  here  we  define  the  more  compact  notation 

TVtW  (n)  :=  v  -  D?  •  (v  -  w)]  n  (3.30) 


to  more  explicitly  indicate  the  dependence  on  the  pre-collision  velocities  and  the  orientation 
angle.  It  should  be  noted  that  the  Jacobian  of  the  transformation  which  maps  (v,  w )  — > 
{V ,w')  is  unity  [30].  Further,  Tv>y  ( n )  =  v.  Babovsky  begins  his  explanation  of  Nanbu’s 
method  by  developing  a  weak  form  of  the  time  discretized  Boltzmann  equation  [8]. 
Rewriting  (3.2)  using  (2.31),  one  obtains 


fk+l  (v)  =  1  -At 


B  (||v  -  w|| ,  6)  f  (w)  dS  in)  dw  fk  (v)  + 


At  f  f  B(||v  -  w\\,6)  f  (V)  fk  (w')dS  (n)dw 

Jr3  Js+ 


(3.31) 


which,  so  long  as  B  is  bounded,  is  an  approximation  to  the  solution  of  (2.28)  with  truncation 
error  on  the  order  of  O  (At). 

Let  (fie  Cb  (R3).  By  multiplying  (3.31)  by  (fi  and  integrating  over  velocity  space,  i.e. 
M3,  one  obtains  a  weak  form  of  the  equation.  Utilizing  the  properties  discussed  above,  it 
can  be  shown  [8]  that  this  expression  can  be  written  as 


where,  Kv<w 


cfi(v)fk+l(v)dv 


Kv,w  ((fi)  f  (v)  f  (w)  dwdv 


(3.32) 


Kvm,  ((fi) 


r 

Js 


1  -  At  B(\\v  -  w\\  ,6)  d6de 


cfi(v)  + 


r 

At  B  (\\v  -  w\\  ,6)  cfi(v')  d6de 

Js 


(3.33) 


is  called  the  transition  kernel.  Note  that  the  kernel  itself  exhibits  no  dependence  on  fk. 
Also,  as  Babovsky  points  out,  there  is  no  guarantee  that  (3.31)  preserves  non-negativity  of 
solutions.  To  ensure  non-negativity  is  preserved  one  must  assume 


1  -  At 


B  (||v  -  w|| ,  6)  dOde  >  0  Vv,  w  e  M3 


(3.34) 
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This  establishes  a  restriction  on  the  time  step  and  typically  also  requires  that  B  be  truncated 
in  velocity  space. 

Babovsky  [8]  proves  the  following  lemma  which  allows  the  transition  kernel  to  be 
written  more  compactly. 

Lemma  3.1.  For  all  v,w  e  M3,  there  exists  a  continuous  mapping,  ®v  w  :  D  — >  S  +,  where 
D  is  the  disk  in  M2  centered  at  the  origin  with  radius  4^,  such  that 

Kv,w  (<P)  =  f  <P (Tv,w  o  cD„  w  (X))  dx  (3.35) 

Jd 

for  any  (p  &  Cb  (M3 ). 

Using  this  property,  and  assuming  B  is  bounded  and  continuous,  we  define  the 
function,  :  D  x  l3  x  I3  ->  K3  as 

Y  (x,  v,  w)  =  TVtW  o  ®,, H,  (x)  (3.36) 

It  can  be  shown  [8]  that  because  B  is  assumed  to  be  bounded  and  continuous,  T  is 
continuous  almost  everywhere  (a.e.)  in  D  x  R3  x  K3.  Physically,  VF  is  a  function  which 
takes  two  velocities  and  a  2-parameter  orientation  as  inputs  and  maps  to  a  single  velocity. 
The  following  definitions  motivate  the  concept  of  a  point  measure  approximation,  which  is 
the  framework  by  which  DSMC  constructs  a  weak  solution  to  the  Boltzmann  equation. 

Definition  3.2.  Let  p  be  a  measure  on  the  measurable  space  (fl,  S).  The  support  of  p  is  the 
set  of  all  points  co  £  Llfor  which  every  open  neighborhood  of  to  has  a  positive  measure  and 
is  denoted  supp  (p ). 

Definition  3.3.  Let  P  be  a  probability  measure  on  the  measurable  space  (Q,  S).  P  is  called 
a  discrete  probability  measure  if  supp  (P)  is  countable. 

Lemma  3.2.  Let  P  be  a  disrete  probability  measure  on  the  measurable  space  (Q,  S)  with 
supp  (P)  =  {cod.  Then  there  exists  {af,  with  a,-  e  M+  such  that  for  any  A  c  fl, 

P(A)  =  Yjai6l0i(A )  (3.37) 

i 
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where  8Wj  is  the  delta  measure  supported  at  <oh  and, 


J]ai=  1  (3.38) 

i 

Definition  3.4.  Let  P  be  a  probability  measure  on  the  measurable  space  (Q,  S),  and  let 
{Pat}  be  a  sequence  of  discrete  probability  measures  on  (Q,S)  with  card  (supp  ( /J  y ))  =  N. 
If  { /Jv }  converges  to  P  in  law,  then  { PN )  is  called  a  point  measure  approximation  to  P. 
Further,  denoting  supp  (YJy)  by  {(oj\p=l,  if 

1  N 

pN(A)  =  -J]s(0i(A)  (3.39) 

i=l 

for  every  A  c  Q,  then  PN  is  said  to  possess  uniform  weights. 

We  also  require  the  following  definitions  and  results  regarding  image  measures. 


Definition  3.5.  Let  (X,  S)  and  (K  IB)  be  measurable  spaces,  and  let  f  be  a  transformation 
from  X  into  Y.  If  f -1  (B)  e  S  for  all  B  e  IB  then  f  is  called  a  measurable  transformation. 

Definition  3.6.  Let  ( X,S,p )  be  a  measure  space  and  (Y.  IB)  a  measurable  space.  Let  T  be 
a  measurable  transformation  from  X  into  Y.  Let  (//  o  T  M  (A)  =  p  (T~l  (A)}  for  all  A  £  IB. 
Since  A  h-»  T~x  (A)  preserves  all  set  operations  and  presents  disjointedness,  p  o  T  l  is  a 
countably  additive  measure  and  is  called  the  image  measure  of  p  by  T.[35] 

Lemma  3.3.  Let  f  be  any  measurable  function  from  Y  into  [-oo,  oo].  Then  ffd  [p  o  j  = 
f  f  o  Tdp  if  either  integral  is  defined,  ( possibly  infinite )  [35]. 

Lemma  3.4.  Let  ( X,S,P )  be  a  measure  space  and  (  Y,  B)  a  measurable  space.  Let  T  be 
a  transformation  from  X  into  Y,  and  let  T  be  continuous  a.e.  Let  j PN  j  be  a  sequence  of 
probability  measures  on  (X,  S)  which  converges  to  P  in  law.  Then,  the  sequence  j  PN  o  T  1  j 
converges  in  law  to  the  image  measure  P  o  T~l  [13]. 


The  following  theorem,  proven  by  Babovsky  [8]  employs  these  results  to  develop 
a  sufficiency  criteria  for  a  point  measure  scheme  to  converge  in  law  to  the  space 
homogeneous  Boltzmann  solution. 
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Theorem  3.1.  Let  P2k  be  the  probability  measure  on  I)  x  M3  x  M3  defined  by 

P2k  ( B )  =  f  f  (v)  f  (w)  dxdwdv  for  any  Be  Dx  M3  x  M3  (3.40) 

Jb 

where  fk  is  the  solution  of  the  space  homogeneous  Boltzmann  equation  at  t  =  4.  Let  Pk+ 1 
be  the  probability  measure  on  K3  defined  by 

Pk+1  (A)  =  I  /+1  (v)  dv  for  any  A  c  M3  (3.41) 

Ja 

Suppose  a  sequence  of  probability  measures  {P2N}  is  a  point  measure  approximation  with 
uniform  weights  to  P2k,  with  support  denoted  by 

supp  (P2n)  =  j(xf ,  vf ,  <)£  (3 -42) 

Then  define  a  probability  measure  PN  on  M3  by 

1  N 

Pn  =  Jj  X  ^(^f.vf.wf )  (A)  /°r  an>’ A  c  r3  (3-43) 

1=1 

where  lF  is  defined  as  in  (3.36).  Then,  {PN}  is  a  point  measure  approximation  to  Pk+l. 


Proof.  By  definition,  card  (supp  (Pn))  =  N.  We  will  show  that  { PN }  converges  to  Pk+1  in 
law.  Take  any  f  e  C/,  (M3).  Then 


('  /f'v  =  .M  l  )  0  (|-  -  'I'  (  >;v.  i".  u,v)j 


(=1 

N 


6  [w  -  wf )  dxdwdv 

f  of(i,  v,  w)  dP2N 


-  Iff 

=  |  f(y)d(P2N  o  T-1)  by  Lemma  3.3 

JlR3  V  7 


(3.44) 
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Thus,  PN  =  P2n  oT^1  .  Recall  that  by  assumption  {P2N}  converges  in  law  to  P2k  as  N  — >  oo. 
Since  is  continuous  a.e.,  by  Lemma  3.4  we  have  {P2N  o  ffi"1!  converges  to  P2k  o  VP_1  in 
law  as  N  — >  oo.  By  (3.44)  this  implies  that  {/V}  converges  to  P2k  o  T*-1  in  law  as  N  — >  oo. 
Next,  for  any  (p  e  Ch  (M3),  we  have 


f  (p{v)d[P2N  of  ')  =  f  f  f  (pUV  (x,v,w))dP2k  by  Lemma  3.3 
J]R3  V  7  Jd  Jr3  Jr3 

m(p  (T  (x,  V,  H'))  f  (v)  fk  (w)  dxdwdv 

3 

nKVtW  UP)  f  (v)  f  (w)  dwdv  by  (3.35)  and  (3.36) 

3 

=  f  cp(v)fk+l(v)dv  by  (3.32) 

JR3 

=  f  (p(v)dPk+i  (3.45) 

JR3 


Thus,  P2k  o>F  1  =  Pk+1  and  hence  we  conclude  that  {PN}  converges  to  Pk+l  in  law  as 
N  -»  oo.  □ 


Theorem  3.1  provides  a  sufficient  criterion  for  the  convergence  of  a  point-measure 
based  simulation  scheme  for  the  space  homogeneous  Boltzmann  equation  in  law.  Namely, 
if  one  can  show  that  a  scheme  constructs  a  point  measure  approximation,  [P2N]  to  P2k  at 
t  =  tk,  then  one  can  construct  a  point  mesure  approximation,  { PN } ,  to  Pk+l  by  defining  its 
support  to  be  given  by  supp  (PN)  =  (supp  ( P2N )). 

Babovsky  [8]  summarizes  Nanbu’s  scheme  more  compactly  than  the  overview 
described  in  the  prior  section.  Babovsky’s  description  lends  itself  better  to  the  proofs 
that  follow  while  retaining  the  important  characteristics  of  the  scheme.  Nanbu’s  scheme 
is  summarized  as  follows: 

•  Let  P{)N  be  a  point  measure  approximation  to  P°  defined  by 

P°  (A)  =  J  f°(v)dv  (3.46) 

where  f°  is  the  initial  condition  for  the  Boltzmann  equation,  and  let  supp  = 
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•  For  a  fixed  N,  choose  a  set  of  N  uniformly  distributed  random  numbers 


rf  G  [0, 1]  and  a  set  of  N  uniformly  distributed  random  vectors  ,  x\  €  D 

Define  collision  partners  C  (. i,N)=  7Vr'v J  +  1  with  velocities  given  by  wf  = 


,N 

C(i,N)' 


Define  the  sequence  of  discrete  probability  measures  P'A1  by 


1  N 

Pn 1  (A)  =  -  J]  <V(*?,vf,w?)  (A)  for  all  Ad3  (3.47) 

i=  1 

This  sequence  will  be  shown  to  be  a  point  measure  approximation  to  P1  defined  by 

P‘(A)  =  Jjl(v)dv  (3.48) 

where  fl  is  the  Boltzmann  solution  at  time  t  =  At 


Repeat  the  process  with  fl  as  the  new  initial  condition. 


Before  continuing,  a  few  additional  definitions  and  results  on  random  variables  and 
probability  theory  are  required. 

Definition  3.7.  If  (fl,  3\,  P)  is  a  probability  space  and  (S /£)  is  any  measurable  space,  a 
function  Xfrom  f 1  into  S  is  called  a  random  variable.  The  law  of  the  random  variable  X, 
denoted  Jl(X),  is  the  image  measure  P  o  X  1  defined  by 

(. P  o  X”1)  (B)  :=  P  (X"1  (B))  for  any  BcS  (3.49) 

The  notation  X  1  (B)  is  used  to  indicate  the  pre-image  of  B,  that  is  X~l  (B)  := 

{co  G  Q. :  X(co)  G  B}.  [35] 


Definition  3.8.  The  expected  value  of  a  random  variable  X  on  (Q,  cr,  P )  is  denoted  EX  and 
is  defined  by  EX  :=  f  XdP.[35] 

Theorem  3.2.  For  any  two  random  variables,  X  and  Y  such  that  EX  and  EY  are  both 
defined  and  finite,  and  any  constant  c,  E  (cX  +  Y)  =  cEX  +  EY. 


57 


Definition  3.9.  The  variance  of  a  random  variable  X  is  denoted  by  var(X )  or  cr 2  ( X ),  and 
is  defined  by 

|  E  (X  -  EX  f  =  EX 2  -  (EX)2  if  EX2  <  oo 
var(X)  :=  a2  (X)  :=  ( 

|  0  otherwise 

If  EX2  <  oo,  then  cr(X)  :=  \J cr1  (X )  is  called  the  standard  deviation  ofX. 

Definition  3.10.  Given  a  random  variable  X  on  a  probability  space  (Q,  S,  P )  with  values 
in  Ma,  (i.e.  X:  Q  — »  ISf),  the  distribution  function  ofX  on  M/f  is  defined  by 


F(x)  :=  P(X  <  x ) 


(3.50) 


where  “<  ”  is  defined  on  Rk  by  x  <y  if  and  only  if  xt  <  ytfor  i  =  1, . . . ,  k. 


These  definitions  motivate  defining  convergence  not  just  for  probability  measures,  but 
for  random  variables  as  well.  In  this  case,  several  forms  of  convergence  are  important. 

Definition  3.11.  Let  { Xn }  be  a  sequence  of  random  variables  with  distribution  functions 
{ F„),  and  X  be  a  random  variable  with  distribution  function,  F.  If  Fn  — >  F,  then  {X,,}  is 
said  to  converge  in  distribution  to  X. 

It  should  be  noted  that  convergence  in  distribution  is  equivalent  to  convergence  in  law 
for  random  variables,  that  is,  Fn  — >  F  if  and  only  if  E  (X„)  —>  E(X)  [12]. 


Theorem  3.3.  (Central  Limit  Theorem  [12])  Suppose  that  {X,,}  is  a  sequence  of 
independent  random  variables  having  the  same  distribution  with  mean  c  and  finite  positive 
variance  cr2.  If  Sn  =  X\  H - +  Xn,  then  the  random  variable 


Sn  -  nc 
cr  sjn 

is  distributed  according  to  the  standard  normal  distribution,  N  (0, 1). 


(3.51) 
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Definition  3.12.  Let  (Q,  S,  P)  be  a  probability  space.  Any  event  with  probability  1  is  said 
to  happen  almost  surely  (a.s.).  A  sequence  of  random  variables  {X,,}  is  said  to  converge 
almost  surely  to  a  random  variable  X  if  P(\oj  g  Q  :  lim,,^  Xn  (to )  =  X  (to)))  =  1.  [35] 

The  following  definition  and  theorem  are  useful  in  considering  the  convergence  of 
stochastic  schemes. 

Definition  3.13.  Given  a  probability  space,  (Q,  S,  P )  and  a  sequence  of  events  {An},  define 
lim  sup  An  as  the  event 

limsupA,,  =  nm>!  U„>„,A„  (3.52) 

The  event  lim  sup  A„  is  also  referred  to  as  A„  infinitely  often  (i.o.),  as  to  G  lim  sup  An  if  and 
only  if  to  G  Anfor  infinitely  many  values  ofn. 

Theorem  3.4.  (Borel-Cantelli  Lemma )  If  {An}  are  any  events  with  P  (A„)  <  oo,  then 
P(limsupA„)  =  0.  If  the  {A,,}  are  independent  and  YmP(An)  =  00  then  P  (limsupA,,)  = 
1  [35]. 

With  these  results,  one  may  proceed  to  prove  the  following  theorem. 

Theorem  3.5.  Let  {PN}  be  a  point-measure  approximation  with  uniform  weights  to  the 
probability  measure  P  on  M3,  defined  by  P  (A)  =  fAf(v)  dv  for  any  A  c  M3.  Further,  let  P 
be  the  probability  measure  on  M3  X  R3  defined  by  P  ( B )  =  fB  f  ( v)  f  (w)  dwdv. 

Denote  supp  (PN)  by  {v,}^.  Define  the  probability  measure,  Pn  on  M3  X  R3  by 

1  n 

PN  (B)  =  —  ^  6VhWi  (B)  for  any  B  gI3x!3  (3.53) 

i=  1 

where  Wi  is  as  defined  in  N anbu’s  method  (i.e.  w,  =  Vcq,N),  where  C(i,N )  =  LXr,J  +  1). 
Then,  PN  is  a  point  measure  approximation  to  P. 

Proof.  Define  the  following  distribution  functions, 

F(v)=  f  f(v')dv'  (3.54) 

Jv'<v 
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G(v,  w)  =  j  j  f  (v')  f  (w')  dv'dw' =  F  (v)  F  ( 

Jv'<v  %Jw'<w 


(3.55) 


which  are  the  distribution  functions  associated  with  the  laws  P  and  P  respectively.  Define 


the  distribution  functions  for  PN  and  PN  by 


r  1  N 

fn(v)  =  I  77  X  5  ( V /  _  v') 

dv'<V  ” 


(3.56) 


F  F  l  \ 

Gn(v,w)  =  I  I  —^^6(v'-Vi)6(w'-Wi)dw'dv'  (3.57) 

*Jv'<v  Jw'<w  z_ j 

respectively.  Recall,  by  Definition  3.11  a  sequence  of  probability  laws  {P,,}  with 
distribution  functions  {F„}  converge  to  a  law  P  with  distribution  function  F  if  and  only 
if  {F„}  converges  to  F.  By  assumption,  we  have  PN  converges  to  P  in  law,  which  implies 
Fn  — >  F.  We  will  prove  that  for  any  ( a ,  b)  e  M3  x  M3,  GN  ( a ,  b)  — »  G  (a,  b )  and  hence  TV 


converges  to  P  in  law.  To  begin,  define  the  following  variables, 

fcw  (v)  =  card  ({z  <  N  :  v,  <  v}) 


(3.58) 


mN(v,w )  =  card  ({z  <  iV  :  (v,-,  wC(/,jv))  <  (v,w)}) 


(3.59) 


Then, 


Pv(v)  = 


kN(v) 


Gn  (v)  - 


N 

mN  (v,  w) 


(3.60) 


(3.61) 


Let  (a,  b)  e 


Then, 


Define 


1  if  v,-  <  a 


0  otherwise 


Gn  (a,  b ) 


mN  (a,  b) 
N 


■  i§  If/'-" 

1  r* 

=  TiX5,  5(w_vc(‘») 

*Jw<b 


8  (v  -  v,)  d(w  -  w,)  JwJv 


(3.62) 
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Let  [i  (s)}^'  be  the  set  of  kN  ( a )  integers  such  that  vl(s)  <  a.  Define 


X, 


Then, 


1  if  vc(i(s)M)  ^  b 


0  otherwise 


kN(a) 


Y  'Vrv 

Gn  ( a ,  b)  =  —  2  Xs 

5=1 


Substituting  the  above  into  (3.61)  we  find 


kN(a) 

mN  {a,  b)  =  ^  Xs 

j=l 


(3.63) 


(3.64) 


(3.65) 


Let  Qn  =  (0, 1).  Then  (X2,„  S„,  A)  is  a  probability  space,  where  S„  is  the  standard  cr-algebra 
on  (0, 1)  and  A  is  the  Lebesgue  measure.  Define  Q  to  be  the  Cartesian  product 


N 


0.  =  Y\ 


(3.66) 


n=  1 


The  elements  of  Q.  are  thus  AMuplcs  {cu„}^=1  with  a>n  e  Q„.  Defining  nm  to  be  the  natural 
projection  of  Q  onto  Q,„  (i.e.  nm  =  a>m),  let  S  be  the  smallest  cr-  algebra  of 

subsets  of  Q  containing  all  sets  n~l  (A)  for  all  n  and  all  A  c  il„.  Here,  n~'  (A)  denotes 
the  pre-image  of  A  (namely,  n~l  (A)  :=  {oj  e  D  :  n(oj)  6  A  c  D„})  .  Then  ( Q,,S,Prob )  is  a 
product  probability  space  with  Prob  =  AN.  Now,  given  an  iV-tuple  {con}f=l  e  Q.  Xs  maps  the 
AMuple  to  either  0  or  1.  That  is,  Xs  :  fl  — >  {0, 1},  which  implies  Xs  is  a  random  variable  on 
the  probability  space  ( Q,S,Prob ).  Furthermore,  the  set  of  random  variables  {Xs}^<ja)  arc 
independent,  and  identically  distributed  (i.i.d.).  Furthermore, 


EXS 


AProb 


fx„ 

f  ...  r  •••  f  xsdAN 

1  •  A  (Q's)  +  0  •  A  (D,  -  Q') 
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where  £2'  =  j.r  e  (0, 1)  :  vL_rWj+i  <  b)  Thus,  NA(Q')  =  card({z  <  N  :  v,  <  /?}).  Therefore, 


EXs  =  T(£2'v) 

na  (n;) 

N 

card({z  <  N  :  Vi  <  b}) 
N 

kN(b) 

N 


Also, 


Var(Xs) 


EX]  -  (EXS)2 


card({z  <  N  :  v,-  <  b})  _  / kN(b)V 
N  \  N  I 

kN(b)  (kNib )\2 

N  \  N  ) 
kN(b)(  kN(b)\ 

N  \  N  ) 


Combining  (3.65)  and  (3.67)  we  have 


kN{a) 

EmN  ( a ,  b)  =  ^  EXS 

s=  1 


kn  ( a )  k„  (b) 
N2 


Define  the  set  (or  event)  AN  c  Q  as 


An  - 


mN  (a,  b ) 

iv 


F  (a,  b) 


>e(N)\ 


(3.67) 


(3.68) 


(3.69) 


(3.70) 
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where  e(N)  is  some  positive  real  number.  In  this  expression,  note  that  the  mN  term  is 
dependant  upon  ,  via  equation  3.65  .  We  next  develop  an  upper  bound  on  Prob  (AN). 
By  the  triangle  inequality, 


(3.71) 


/( r.\N  . 

mN  ( a ,  (?)  Fv  (a)  kN  (b)  ^ 

N  N  N 

2  I 

|  trdi=i  • 

((Of.,  : 

kN(a)kN(b ) 

iV  Af 

2  i 

Thus, 


Prob  (An)  <  Prob\^[ri}i=l  : 

Probii^}^  : 


mN  ( a ,  b)  kN  ( a )  kN  (b)^  ^  e  (N) 

N  N  A T*  >  2 


I  kN  (a)  kN  ( b ) 


-F(a)F(b)\  > 


(3.72) 


U'  ",=1  |  N  N  \  2  }) 

Notice  that  by  equations  3.56  and  3.60,  the  term  _  F  (a)  F  (b) |  does  not  depend 

on  and  is  simply  a  constant  for  a  given  N.  Define  e2  (N)  to  be 

_  ^  \kN  (a)  kN  (b)  r,_Anm| 


e2(N)  =  2 


N  N 


F  (a)  F  (b)\ 


and  thus, 


Prob  \  {r,(.=1  : 


I  kN  {a)  kN  (b) 


F  (a)  F  (b)\  > 


e2(N) 


(3.73) 


(3.74) 


((/  ",=  '  |  N  N  |  2  }) 

Note  that  since  PN  — >  P  implies  FN  — >  F,  we  have  e2  (/V)  — >  0  as  iV  — >  oo. 

Since  mN  is  the  summation  of  the  kN  i.i.d  random  variables  Xs,  by  the  central  limit 

theorem  (Theorem  3.3),  we  have  that 

mN  ( a ,  b )  -  kN  ( a )  EXS 
VFv  ( a )  VarXs 

is  distributed  according  to  the  standard  normal  distribution  N  (0, 1).  This  implies  that 


-  £  (X,  -  EXs)  = 

5=1 

is  distributed  according  to  N  (0,  crN),  where 

VFv 

cr  N  -  - 


mN  (a,  b)  kN  ( a )  kN  (b) 

N  N  N~ 


(3.76) 


/c,v  (a)  Var  (Xv) 

iv 

lkN(a)  kN(b)  ( kN(b)\ 

N  N  X1  N  ) 

Viv 


(3.77) 
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Now,  since  E 


m^(a,b)  k^{a)  k^(b) 

N  ~N  N~ 


j  =  0,  and  crN  is  finite,  choose  e,  (N)  to  be  the  smallest 


positive  number  such  that 


n  mN(a,b)  kN(a)kN(b )  e,  (N) 

U  N  N  N  2 


(3.78) 


The  existence  of  such  a  value  follows  from  the  conditions  on  the  expected  value  and 
variance  of  mN(“’h)  -  Furthermore,  ej  (N)  — >  0  as  N  — >  oo  since 


lim  crN  =  lim 

N—>oo  N—>co 


j  kN(a )  kN{b) 
N  N 


^(i-T) 

Viv 


=  lim  —  lim  A  — 

\V^oo  J  I  /V-> OO  V  N  N 


kN  (a)  kN  (b)  /  _  kN(b) 
N  N  \  N 


(0)  JF(a)F(b)(l-F(b)) 


Let  e  (AO  =  max  {ei  (N) ,  e2  (AO).  Then,  e  (N)  — >  0  as  N  co  and 


mN  (a,  b )  kN  ( a )  kN  (b)  e  (N) 
N  N  N~  >  2 


Prob  (An)  <  Prob  ( {r,}J=1  : 


\  T7  1 1  \  ^  fW 


+Prob  ^|{r;}7  j  ;  -  -  ^  -  F  (a)  F  (b)  > 

=  Probl(lr)N  •  k^a)kN(b)  e(N) V\ 

l|l,,,=1'  N  N  N  2| 


Thus,  we  have, 


OO  OO  J 

YjProb{AN)-YjW<0° 


(3.79) 


Thus  defining  A  to  be  the  event  AN  i.o.,  as  given  in  (3.52),  by  the  Borel-Cantelli  Lemma 
(Theorem  3.4)  we  have  Prob  (A)  =  0.  Thus, 


Urn  m"(a-b)-F(tt)Fib)  =  o 

oo  N  v  v 


(3.80) 


for  almost  all  r.  Combining  this  result  with  (3.55)  and  (3.61)  yeilds  that 


lim  \Gn  (a,  b)  -  G  (a,  b)\  -  0  almost  surely 

N — ^oo 


(3.81) 
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for  all  ( a ,  b)  e  M3  x  M3.  Hence  PN  converges  in  law  to  P  as  N  — »  oo  almost  surely.  □ 

The  next  result  incorporates  the  distribution  of  collision  orientations;  the  proof  will  be 
omitted  as  it  follows  the  same  technique  as  Theorem  3.5. 

Theorem  3.6.  Let  { /Jy }  be  a  point  measure  approximation  to  the  probability  measure  P  on 
7i  \  defined  by  P  (A)  =  fAf(v)  dv  for  any  A  e  7d\  Further,  let  P  be  the  probability  measure 
on  Dxt3xR3  defined  by  P  ( B )  =  JBf(v)f  ( w )  dxdwdv  for  any  B  e  D  xl3  x  M3. 

Denote  supp(Pw)  by  {(v,,  Let  be  a  sequence  of  uniformly  distributed 

random  numbers  on  D  =  j.r  6  M2  :  ||v||  <  -4=  j,  where  ||-||  represents  the  Euclidean  norm. 
Define  the  probability  measure  Pn  on  D  X  M3  X  M3  by 

1  n 

PN  (B)  =  —  SXi,Vi,wi  (B)  for  any  6eDxR3xI3  (3.82) 

1=1 

Then  PN  is  a  point  measure  approximation  to  P. 

Theorem  3.7.  Nanbu ’s  scheme  converges  in  law  almost  surely  to  the  solution  of  the  time 
discretized  space  homogeneous  Boltzmann  equation  (3.2)  for  all  time  steps. 

Proof.  Combing  Theorems  3. 1,3-5,  and  3.6  yields  the  desired  result.  □ 

3.4  Nanbu’s  DSMC  Scheme  Applied  to  the  Bobylev-Krook-Wu  Problem 

To  provide  a  numerical  baseline  for  comparison  with  the  schemes  developed  by  the 
author,  Nanbu’s  DSMC  technique  was  applied  to  the  Bobylev-Krook-Wu  Problem.  As 
Nanbu’s  scheme  converges  only  in  law,  it  is  not  possible  to  directly  evaluate  the  error 
present  in  the  solution  (e.g.  in  terms  of  L1  (M3)  norm,  etc.).  For  this  reason,  the  moment 
expressions  for  the  Bobylev-Krook-Wu  solution  presented  in  Chapter  2  provide  a  very 
valuable  and  commonly  exploited  tool  for  evaluating  the  effectiveness  of  various  DSMC 
schemes.  Utilizing  the  Nanbu  approximation  to  the  distribution  function,  these  normalized 
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moments  can  easily  be  computed  as 


Zn  (t) 


l  Np  i  r 

(2/t  +  l)!!  §  iVp  L 


d(v-Vim\v\\2ndv 


1 


N„ 


Np  (2n  +  1) 


M 


(3.83) 


'P'-~  ■  i=i 

Figure  3.1  presents  a  comparison  of  the  first  four  normalized  moments  of  the  density 
function  for  the  case  in  which  Np  =  100.  The  results  have  been  averaged  over  an  ensemble 
of  600  runs.  This  is  common  practice  to  reduce  the  variance  present  in  the  solution. 
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Figure  3.1:  Normalized  Moments  of  Bobylev-Krook-Wu  Solution  (Np  =  100,  600 
run  ensemble) 

Notice  the  significant  amount  of  variance  present  in  the  solution  and  that  the  solution  tends 
to  diverge  more  significantly  for  the  higher  order  moments.  A  particularly  concerning 
effect  for  which  Nanbu’s  method  has  received  much  criticism  is  apparent  in  examining  z,\ . 
Note  that  the  exact  solution  for  z\  is  unity  for  all  time.  This  moment  is  in  fact  a  multiple 
of  the  energy  in  the  gas  which  must  be  conserved  for  all  time.  Nanbu’s  method  does  not 
exactly  conserve  energy  with  each  collision,  but  rather  only  on  average.  This  leads  to  some 
variation  in  z\  throughout  the  simulation. 

The  results  illustrated  in  Figures  3.2  and  3.3  illustrate  results  for  the  first  and  second 
normalized  moments  for  varying  values  of  Np.  These  results  were  derived  from  an 
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ensemble  of  400  samples,  and  the  time  step  was  chosen  as  At  =  l/Np.  Notice  that  as  the 
number  of  simulated  particles  is  increased,  both  the  accuracy  and  variance  are  improved; 
however,  even  at  the  largest  number  of  particles  simulated,  the  variance  is  still  substantial. 
Higher  order  moments  exhibit  similar  features,  however  the  divergence  from  the  exact 
solution  increases  as  higher  moments  are  considered. 


t 

Figure  3.2:  First  Normalized  Moment  with  Varying  Np  (400  run  sample) 


To  measure  the  effect  of  the  simulation  parameters  on  the  variance  in  the  solution,  one 
technique  is  to  consider  the  total  variation  of  each  moment.  For  the  current  problem  we 
define  the  total  variation  of  the  nth  moment  as 

rtfinal 

V  (z.n)  =  \Zn(t)\dt  (3.84) 

Jo 

For  the  current  work,  tfinaj  =  10.  The  total  variation  for  Np  =  100  is  plotted  against  the 
number  of  runs  in  the  ensemble  average  in  Figure  3.4  along  with  the  exact  solution  for 
the  total  variation  of  each  moment.  Notice  that  the  only  moment  whose  total  variation  is 

identically  zero  for  all  time  is  Zi  which  again  is  a  statement  of  the  conservation  of  energy. 
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Figure  3.3:  Second  Normalized  Moment  with  Varying  Np  (400  run  sample) 


Figure  3.4:  Total  Variation  of  Nanbu  Method  as  a  Function  of  Sample  Size  for 
Np  =  100 


While  increasing  the  number  of  samples  in  a  DSMC  simulation  is  the  primary 
approach  for  variance  reduction,  it  does  not  improve  the  accuracy  of  the  solution  itself.  To 
improve  the  accuracy  of  a  simulation  one  must  increase  the  number  of  simulated  particles. 
Since  it  is  not  possible  to  directly  analyze  the  distribution  function,  consider  the  L1  (lR3J 
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error  for  each  of  the  moments  defined  by 


Error  (zn) 


S'* final 

=  1 


I Zn  (t)  ~  Zn  (01  dt 


(3.85) 


Fixing  the  number  of  samples  at  1000,  Figure  3.5  shows  the  effect  of  increasing  the  number 
of  particles  on  the  L1  (iR3)  error  in  the  solution. 


Figure  3.5:  L1  (lR3J  Error  in  Normalized  Moments  as  a  Function  of  the  Number  of 
Simulated  Particles,  Nanbu  Method  (1000  run  ensemble) 
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IV.  A  Distributional  Monte  Carlo  Algorithm  Based  on  Kernel  Density  Estimation 


As  a  first  attempt  towards  a  distributional  approach,  the  concept  of  Kernel  Density 
Estimation  was  applied  to  the  DSMC  framework.  As  outlined  in  Chapter  1,  the 
Distributional  Monte  Carlo  method  advocated  by  the  author  replaces  the  point  measure 
representation  of  the  velocity  density  function  with  one  in  which  each  simulated  particle  is 
allowed  to  possess  an  entire  velocity  density  function  versus  only  a  single  velocity  vector. 
One  interesting  limiting  case  is  the  case  in  which  particle  velocity  density  functions  are 
all  assumed  to  have  the  same  functional  form.  Perhaps  the  simplest  such  representation  is 
obtained  by  writing  the  particle  density  functions  as 

m  =  kK[nr)  (4'1) 

where  A-  is  a  specified  probability  density  function,  h  e  M+,  and  v,  is  the  mean  velocity 
of  the  ith  simulated  particle.  The  problem  then  becomes  one  of  determining  a  suitable 
stochastic  scheme  for  altering  v,  to  account  for  the  collisional  process.  In  this  chapter, 
such  a  scheme  is  developed,  which  is  equivalent  to  applying  kernel  density  estimation  to 
Nanbu’s  DSMC  scheme.  It  will  be  shown  that  with  an  appropriately  chosen  value  of  K , 
the  scheme  converges  in  law,  as  well  as  in  solution  for  L°°  (iR3)  and  bounded  solutions  of 
the  space  homogeneous  Boltzmann  equation.  These  latter  two  forms  of  convergence  have 
never  previously  been  demonstrated  for  a  stochastic  particle  scheme.  The  method  could 
therefore  be  viewed  as  a  bridge  from  Direct  Simulation  to  Direct  Solution. 

The  remainder  of  this  chapter  is  outlined  as  follows.  First,  a  brief  explanation  of  kernel 
density  estimation  is  presented,  as  the  development  of  the  scheme  draws  significantly  upon 
the  concept.  Next,  a  brief  derivation  of  the  technique  is  presented,  followed  by  proof  of  its 
convergence  in  law,  as  well  as  in  solution  for  L°°  (iR3)  and  bounded  solutions  of  the  space 
homogeneous  Boltzmann  equation.  The  chapter  concludes  with  some  numerical  results 
obtained  by  applying  the  technique  to  the  Bobylev-Krook-Wu  solution. 
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4.1  Overview  of  Kernel  Density  Estimation 


Kernel  Density  Estimation  (KDE)  is  a  technique  for  estimating  the  probability  density 
function  of  a  random  variable  X  e  W1  from  a  set  of  discrete  samples  as  follows  [77]. 

'x-X, 


N 


f  (x;  h) 


Ly 

hd  4_J 

i—  1 


Nhd 


K 


h 


(4.2) 


Here,  N  is  the  number  of  discrete  samples,  X,  the  value  of  the  ith  sample,  he  M+  the 
kernel  bandwidth,  and  Ke  L2  the  kernel  function.  The  kernel  function  must  satisfy 
the  following  conditions 


K  ( x ) dx  =  1 

d 

xK  (x)  dx  =  0 


(4.3) 

(4.4) 


The  problem  then  becomes  one  of  determining  a  suitable  K  and  h  with  which  to 
approximate  the  distribution  function.  The  value  of  h  is  chosen  to  minimize  the  error 
between  the  estimator  and  the  actual  distribution  function  in  some  sense.  If  h  is  too  small, 
the  estimator  will  exhibit  overly  oscillatory  behavior.  It  h  is  too  large,  subtle  features  of 
the  distribution  function  may  not  be  captured  by  the  estimator.  Wand  [77]  shows  that  the 
asymptotic  mean  square  error  between  /  and  /  is  minimized  when 


m  (K) 

M  ( K))2m(f")N_ 


(4.5) 


where, 


m  ( g ) 


d2  ( 8 ) 


=  f  [g(x)f 

jRd 

f 


dx 


x2 g  (x)  dx 


Notice  that  calculation  of  such  an  h  requires  not  only  that  f"  is  known,  but  also  that 
/  6  W2-2.  If  /  is  normal  with  variance  cr2,  (4.5)  becomes, 

8  \fnin  (K) 


h 


3<ju2(K))2N\ 


cr 


(4.6) 
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4.2  Application  to  Distributional  Monte  Carlo  (DMC-KDE) 

Observe  that  a  point  measure  approximation  (as  in  DSMC,  e.g.  (1.11))  to  the 
distribution  function  may  be  viewed  as  a  special  case  of  a  kernel  density  estimator  with 
K  =  8  and  h  =  1.  Recognizing  this  similarity,  a  distributional  Monte  Carlo  method 
employing  some  of  the  results  from  kernel  density  estimation  was  developed  by  the 
author  [66,  67,  69].  The  approach  has  been  termed  DMC-KDE.  It  should  be  noted  that 
KDE  has  been  applied  in  the  Variance  Reduced  DSMC  (VRDSMC)  approach  proposed  by 
Al-Mohssen  and  Hadjiconstantinou  [3-5]  though  in  that  application  it  was  employed  only 
as  a  smoothing  and  stabilization  technique. 

In  the  Distributional  Monte  Carlo  approach  we  allow  each  simulated  particle  to 
possess  its  own  velocity  distribution  function,  f.  The  overall  distribution  function  in  the 
gas  is  then  given  by 


In  the  DMC-KDE  approach,  we  make  the  simplification  that  each  particle’s  velocity  is 
distributed  according  to  a  prescribed  distribution,  that  is 

m  =  vK(nr)  <4-8) 


The  mean  of  each  particle’s  distribution  function  v,  is  allowed  to  vary,  but  the  kernel 
function  and  bandwidth  are  chosen  to  be  identical  for  all  particles.  Therefore  the 
approximation  to  the  overall  distribution  function  of  the  gas  is  given  by 


npn  i= i 


v  -  Vi 

h 


(4.9) 


Choosing  h  :  N 


with  the  property 


lim  h  (np)  =  0  (4.10) 

Np  — ^oo  V  ’ 

and  choosing  K  with  properties  as  described  in  the  previous  section,  (4.9)  is  observed  to  be 
a  kernel  density  estimator  for  /.  Although  we  prove  that  convergence  is  guaranteed  for  any 
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such  h  and  K,  it  is  necessary  to  choose  specific  values  of  such  parameters  from  which  to 
construct  a  simulation  scheme.  On  the  basis  of  physical  reasoning  we  choose  a  Gaussian 
kernel  for  K  and  utilize  (4.6)  to  determine  h. 

K(x)  =  (2;r)-3/2exp(-^J  (4.11) 

1 

5 

crest  (4.12) 

where  crest  is  an  estimation  of  the  standard  deviation  of  /.  These  choices  are  advantageous 
for  a  number  of  reasons.  Since  h  is  chosen  such  that  h  —>  0  as  Np  — >  oo,  (4.9)  will 
converge  to  the  delta  representation  when  Np  becomes  large,  recovering  the  point  measure 
approximation  of  DSMC.  Further,  the  distribution  function  of  each  simulated  particle  is 
Maxwellian,  the  prevailing  distribution  in  an  equilibrium  gas.  The  physical  interpretation 
therefore  is  that  although  the  collection  of  particles,  which  a  simulated  particle  represents, 
all  possess  different  velocities,  as  a  collection,  the  particles  represented  by  a  single 
simulated  particle  are  in  translational  equilibrium  with  one  another.  This  represents  a 
relaxation  of  the  assumption  made  by  DSMC  that  the  collection  of  particles  possess  the 
same  singular  velocity.  To  develop  the  mathematical  formulation  of  this  approach,  we 
follow  an  analysis  similar  to  Nanbu  [58]. 

Beginning  with  (4.9)  we  seek  to  determine  the  evolution  of  the  distribution  function 
due  to  intermolecular  collisions  through  the  time  interval  At.  We  begin  by  utilizing  a 
forward  Euler  discretization 


\  _ 

32 

)  - 

3V2 Np 

f(v,t  +  At)  =  f  (v,  t)  +  A t^f  (v,  t) 

at 

where,  'jj-  is  obtained  from  the  space  homogeneous  Boltzmann  equation 


%(v,t)  =  Q(f,f)(v,t) 
ot 


(4.13) 


(4.14) 


and  Q  is  defined  by  (2.31).  Substituting  (4.9)  into  (4.14),  one  obtains 


It  ~  T‘>) 


(4.15) 


P  «=1  7=1 
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where, 


Sij 


=*Li 

'4Ii 


,V  -  V;  w  ~  Vi  , 

K\^^\K\ — — |  •  \\v-w\\B(\\v-w\\,6)dndw 


k['^l)k  ‘w-v’ 


h  I  \  h 

Substituting  (4.11),  (1.7)  and  (1.8)  into  (4.16),  one  obtains 

=  2^lG(V’W)eXp 


-  nil  B  (||v  -  nil ,  6)  dndw 


Sij 


k  [li-"2 '  “""2 


4  h2 


+  v  r+ 


+ 


|v,||-  -  (v,  +  V,)  •  (w  +  v)  |  ||v  -  vv 1 1  dw 


where, 


G  (v,  w) 


L 


exp  [a  ■  n]  B(]\v  -  w\\  ,9) 
\\v-w\\  (Vj  -  V/J 


4  h2 


Nanbu  [58]  shows  that  G  may  be  approximated  for  small  h  by  the  following 

eM 

G  (v,  w)  ~  2 nB  (||v  -  nil  ,x)  —  for  h  small 

Hall 

Substituting  this  expression  for  G  into  (4.18)  yields 


Stj  = 


|v;-  -  Vj\ 


/ 

Jm3 


f  (2/r) 


■3/2 


/l3 


exp 


I V  +  W  -  Vj  -  Vj  j 


\  1 

(||v  —  w|  |  —  gjjj" 

p 

Ah 2 

4/i2 

B(\\v -w\\,x)  dw 


(4.16) 

(4.17) 


(4.18) 


(4.19) 


where  gij  =  ||vy-  -  v,||.  Next,  we  consider  (4.13)  under  the  limit  as  h  — >  0  with  ^  as  in 
(4.15). 

.  Np  N„ 

f  (v,  t  +  At)  =  f  (v,  ,)  +  t2  £  £  (s i  -  Tj)  (4.20) 

>  i=  1  j=  1 

where, 


S7,  =  limS,, 

/i->0  ; 

r*.  =  lim  Tij 

J  /!-» 0  J 


1A 


Definition  4.1.  Let  {#„}  be  a  family  of  locally  integrable  functions  on  R"  with  parameter 
if  el.  {#„}  is  called  an  /(-dimensional  delta  family  as  a  — >  a0  if 

lim  I  ga  (a:)  (p  O)  dx  =  (p  (0) 

a^a°  JR" 

where  (p  is  any  bounded  continuous  function  on  W.  We  write 


lim  ga  (*)  =  5  (.r) 

a— >ao 


in  conformance  with  [72]. 


For  any  bounded  and  continuous  function  (p  on  M3,  it  can  be  shown  that  It  can  be 
shown  that 


C  I  (2n)~3/1 

tnjJ,R<‘R_^exp 


(v  +  W -Vi  -  Vj} 


4h2 


dw  =  (p  (vj  +  Vj  -  v)  (4.21) 


Thus,  the  first  bracketed  term  in  (4.19)  is  a  three-dimensional  delta  family.  Utilizing  this 
property,  it  can  be  shown 


s*ij  =  ~5(l|yt|1  _  l'8ij}<r(sij>x) 


where, 


v  =  v  ■ 


^(u-vy) 


Using  (4.17)  and  performing  a  similar  analysis,  it  can  be  shown 


(4.22) 


(4.23) 


T*j  =  gijBijS  (v  -  Vj) 


(4.24) 


where, 


Bi 


dn 


Substituting  these  terms  into  (4.20),  one  obtains 


1  Np 

f  (v,  t  +  At)  =  —  V  [(1  -  Pi)  6  (J-  \’i)  +  Qi  (v)] 

np  i=i 


(4.25) 
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where, 


At  ^ 

P'  ~  Tv  X  gijBij 


4Ar  y  B{gii,x) 
N  ^  2- 


^lllvll-  L,7 


(4.26) 

(4.27) 


and  x  is  the  angle  between  v*  and  v,  -  vr  Pl  represents  the  probability  that  the  ith  particle 
collides  in  the  time  interval  At,  while  the  individual  terms 


Pij  ~  pj8ijBU 

represent  the  probability  that  the  ith  particle  collides  with  the  jth  particle  over  At.  Also, 
notice  that  Qt  represents  the  portion  of  distribution  function  describing  the  effects  of 
collisions  over  the  time  interval  At.  Having  passed  to  the  limit  of  large  Np  (where  the 
distributions  tend  towards  a  delta  approximation),  we  have  obtained  the  same  result  as 
Nanbu  [58]  for  the  evolution  of  the  distribution  function  over  At.  Therefore,  we  may 
reuse  the  collision  selection  and  modeling  rules  developed  by  Nanbu,  but  with  a  new 
interpretation.  Namely,  v,  now  represents  the  mean  velocity  of  the  ith  simulated  particle. 
Collision  interactions  therefore  have  the  effect  of  shifting  the  individual  Maxwellian 
densities  to  new  mean  values.  The  stochastic  scheme  to  evolve  /  through  At  is  therefore 
given  as  follows: 


•  For  each  particle,  calculate  P(.  Generate  a  random  number  in  the  interval  (0, 1) .  If 
Pi  >  r\,  accept  the  particle  for  collision. 

•  Sample  a  collision  partner  j  from  the  conditional  probability  distribution  P*jk  = 

by  sampling  a  second  random  number  r2  uniformly  from  the  interval  (0, 1)  and 
identifying  the  j  which  satisfies  P*k  <  r2  <  Y,i=  \  P*k - 

•  Sample  the  direction  of  v*  based  on  (4.27),  and  compute  the  post  collision  velocity 
of  the  ith  particle. 
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Like  DSMC,  the  simulation  would  be  evolved  many  times  to  generate  an  ensemble 
averaged  solution  so  as  to  reduce  statistical  fluctuations.  Nanbu’s  scheme  has  in  the  past 
been  criticized  for  not  mainintaining  strict  conservation  of  energy  in  each  collision  but 
only  over  the  ensemble.  Since  the  current  scheme  employs  a  similar  sampling  procedure 
for  the  Maxwellian  centers,  it  will  not  conserve  energy  with  each  collision  either.  This  is 
not  a  major  concern  in  demonstrating  the  benefits  of  such  an  approach,  and  a  number  of 
practical  alternative  sampling  procedures  for  collision  interactions  exist  which  do  conserve 
energy  with  each  collision. 

Although  the  scheme  is  in  some  sense  similar  to  Nanbu,  the  effect  of  allowing 
velocities  to  be  distributed  has  a  significant  impact  on  the  mathematical  convergence 
properties  of  the  method.  Namely,  whereas  DSMC  can  only  achieve  convergence  in 
probability  measure  (weak  convergence),  it  is  proven  that  the  DMC-KDE  approach  results 
in  convergence  in  solution  (strong  convergence)  for  L°°  (R3J  and  bounded  solutions  of 
the  Boltzmann  equation.  Therefore,  rather  than  a  stochastic  simulator  of  the  Boltzmann 
equation,  the  DMC-KDE  approach  represents  a  stochastic  solver  of  the  Boltzmann 
equation. 


4.3  Proof  of  Convergence  of  DMC-KDE  Approach 

In  the  following,  we  prove  weak  convergence  of  the  DMC-KDE  approximation  for 
the  space  homogeneous  Boltzmann  equation.  Based  upon  the  results  of  Babovsky  and 
Illner  [8,  9]  this  is  not  unreasonable  to  expect.  Although  in  the  previous  section  specific 
functions  for  h  and  K  were  chosen,  the  proof  is  for  the  more  general  case. 

_  N 

Let  {v,)(=j  be  the  mean  velocities  of  the  Np  simulated  particles  at  a  given  time  step 
derived  by  the  above  method.  The  velocity  distribution  function  (VDF)  of  the  DMC-KDE 
method  is  then 


/(v)  = 


1  Np 

tt^k 


V  -  Vj\ 

h  I 
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where 


geL2(l3):g(r)>0Vr;  f  g(x)dx=  1; 

{  Jr3 

Xxg  (x)  dx  -  0  > , 
3  J 

and  h  :  N  — »  M+  defined  as  in  (4.10).  Define  {//,}  by 


Then,  /  may  be  rewritten  as, 


Mx) = lAl 


1 

/(v)  =  ^  fh  (v  -  V,-) 


p  i=l 

Lemma  4.1.  {//,}  is  a  delta  family  as  h  — >  0+. 


Proof.  [72]  Let  u  =  f.  we  have 


f  fh  (x)dx  =  —  f  K(-)dx  =  f  K  (u)du  =  1 . 
Jr3  «  Jr3  V«/  Jr 


Also,  for  any  A  >  0, 


lim  f 

,,_>0  J||*||>A 


lim  fh  (x)dx  =  lim 


imif 

,_>0  J\\x\\>A 


K  - )  dx 


lim  I  K(u)du 

h^°  J\\u\\>i 

0 


and, 


h->  0 


/  ■ 

J||x||<A 


lim  I  fh  (x)  dx  =  lim 


lix i 


^  I  -Wr 


*-»°  ^3  Jimka  “  V* 

=  lim  I  K(u)du 

h^°  J||«ll<f 

=  1 


Now,  let  f  be  any  bounded  and  continuous  function  on  M3.  We  have, 


lim  f 

h~ *°  Jr 


lim  fh  (x)  (p  (x)  dx  -  f  (0)  =  lim 


fh  (x)  f  (x)  dx  -cp{ 0)  fh  (x)  dx 


f 

Jr3 

lim  f  fh  (x)  [(p  (x)  -  (p  (0)]  dx 
h~*Q  Jr3 


f 

Jr3 


(4.28) 
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Define  rj(x)  =  f  (x)  -  f  (0).  Let  e  >  0  and  B  >  0.  We  have, 


fh  (x)  7]  (x)  dx 


I  fh  (x)  n  (x)  dx  + 
Jm\<b 


I  f,  (x)  rj  (x)  dx 
J\\x\\>B 


Choose  M  >  0  such  that  f(x)\  <  M  V.r.  Let  p  (B)  =  max||x||<B  \rj(x)\.  We  have, 


/ 

Jr3 


fh  (x)  rj  (x)  dx 


< 


I  fh  (x)  ij  (x)  dx 
J\\x\\<B 

I  fh  ( x )  dx 

J\\x\\<B 


<  P  ( B ) 

<  p  (B)  +  M 


+ 


+  M 


I  fh  (x)  T]  (x)  dx 

J\\x\\ >B 

I  fi,  (x)  dx 
J\\x\\>B 


I  fh  (x)  dx 
J\\x\\>B 


Since  r/  is  continuous  and  rj  (0)  =  0,  there  exists  6eR  such  that  p  ( B )  <  Also,  from 
above  we  have  that  there  exists  a  >  0  such  that 


/  . 

J\\x\\>B 


fh  (x)  dx 


< 


2  M 


whenever  0  <  h  <  a.  Therefore,  for  any  e  >  0,  we  have  shown  that  there  exists  a  such  that 


f 


fh  (x)  rj  (x)  dx 


<  p  (B)  +  M 


I  fh  (x)  dx 
J\\x\\>B 


<  -+M— 
2  2  M 


whenever  0  <  h  <  a.  Therefore,  lini/^o  fh  (x)  rj  (x)  dx  -  0.  Which  implies 
lim/,^o  fh  (x)  (f)  (v)  dx  =  cp  (0),  and  hence  {//,}  is  a  delta  family  as  h  —>  0.  □ 

Recall  that  the  parameter  h  is  chosen  in  the  DMC-KDE  method  such  that  limiV()^co  h  ^Np J 
0.  Then  {//,(a?;))}  is  a  delta  family  as  Np  — >  oo  by  Lemma  4.1.  Denote  this  family  by  {  /,!V(J. 
Next  it  is  proven  that  in  the  limit  as  Np  — >  oo  the  probability  measure  generated  by  the 
DMC-KDE  approach  is  the  same  as  that  generated  by  the  Nanbu  DSMC  method. 
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Lemma  4.2.  For  any  bounded  and  continuous  <fi  on  M3, 

lim  (  <p(v)  f  (v)  dv  =  lim  (  f  (v)  /  (v)  dv 

J R3  Np^°°  J R3 

where, 

1 

/(v)  =)—  V  5 (v  -  Vi)  (4.29) 

,=i 


Proof.  Choose  any  e  >  0.  Then  by  Lemma  4.1,  for  any  bounded  and  continuous  f,  there 
exists  M  such  that 


0  (v)  /n„  (v  -  V,-)  dv  -  f  (v,-)  <  e 

Ur3 

for  all  >  M .  Recalling  that  the  DMC-KDE  scheme  generates  the  same  values  for  the 
Maxwellian  centers  as  the  Nanbu  method  generates  for  molecular  velocities,  we  have 


f  0(v)(/C 

Jr3 


v)  -  /  (v 


XI  X  ^ 

0(v)—  X  /vp  (V  -  V/)  -  <S  (v  -  V/)  dv 

3  'v  Itr  J 

<  X  J"  0  (v)  (/vy  (v  -  V/)  -  <5  (v  -  V/))  dv 

<  s:2> 


for  all  Np>  M.  □ 

Combining  this  result  with  Theorem  3.7,  yields  the  following  result. 

Theorem  4.1.  If  the  time  discretized  space  homogeneous  Boltzmann  equation  (3.2)  has 
a  non-negative  solution  f  e  L1  (M3J,  then  the  solution  f  4.9  of  the  DMC-KDE  method 
converges  in  law  (see  Definition  3.1)  at  each  time  step. 

Proof.  Take  any  f  e  C/,  (K3).  Choose  any  e  >  0.  Then,  by  the  triangle  inequality, 

f  0(v)(/(v)  -f(v))dv  <  f  <p(v)(f(v)  -f(v))  dv 

Jr3  Jr3 

+  f  <P(v)(f(v)  -/(v))  dv 

Jr3 
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where  /  is  as  defined  in  (4.29).  Applying  Lemma  4.2  and  Theorem  3.7  to  the  terms  to  the 
right  of  the  inequality  yields  the  desired  result.  □ 


We  have  therefore  proven  that  the  DMC-KDE  method  exhibits  the  same  convergence 
as  Nanbu’s  method  in  the  general  case,  namely  convergence  in  law  to  the  space 
homogeneous  Boltzmann  solution.  We  next  prove  that  stronger  forms  of  convergence  are 
possible  compared  with  Nanbu’s  method,  specifically  for  solutions  which  are  L°°  (R3J  or 
bounded.  Such  solutions  arise  frequently  in  kinetic  theory  and  are  of  greater  practical 
interest  than  L1  (M3)  solutions. 

Corollary  4.1.  If  the  time  discretized  space  homogeneous  Boltzmann  equation  (3.2)  has 
a  non-negative  solution  f  6  L°°  (R3),  then  the  solution  f  (4.9)  of  the  DMC-KDE  method 
converges  in  L°°  (R3J  to  f  at  each  time  step.  That  is, 


lim  11/ -/II  =0 

N„->  CO  11  1,00 


Proof.  Take  any  e  >  0.  Since  /,  /  e  L°°  (R3),  there  exist  B \ ,  B2  s  BA  such  that  |/(jc)|  <  //, 
and  |  /( x)\<  B2  almost  everywhere.  Let  S  i  and  Si  be  the  sets  of  zero  measure  over  which 
these  inequalities  do  not  hold  for  /  and  /  respectively.  Let  S  =  S  |  U,S'2.  For  any  x'eR3A, 
define  :  M3  — >  M  by 


fh  U) 


1 


exp  - 


\x  -  x 


,/||2 


s/2nh3 

where  h  £  RA.  By  Lemma  4.1,  cf)h  is  a  delta  family  as  h  0 
there  exists  H\  >  0  such  that 


centered  at  x' .  Therefore, 


r 

JR3 


fh(x)f(x)dx-f(x') 
for  all  h  <  Hi.  Likewise,  there  exists  H2  >  0  such  that 

f  <Ph(x)f(x)dx-f(x') 

J  R3 


e 

<  - 
3 


6 

<  - 

3 


(4.30) 


(4.31) 
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Note  that  (ph  is  everywhere  continuous  and  \<ph  (x)|  <  ^=- h3  for  all  x  €  M3 .  Thus  by 
Theorem  4.1,  there  exist  N  >  0  such  that 


6 

<  - 
3 


J'  (ph  O)  (I(x)  -  f  (x))  dx 
for  all  Np  >  N.  Choose  H  <  min  {Hi,  H2),  then  for  any  Np  <  N  we  have, 


(4.32) 


\\f~f\L  =  esssup|/-/| 

=  sup  \f(x')-f(x')\ 

fix')-  f  (phix)fix)dx 

J  R3 


<  sup 

x'sK  3-S 


+ 


r 

Jr3 


(phix)fix)  dx-  fix) 


<  sup 

x'sR  3-S 


fix')-  f 

Jr3 


by  the  triangle  inequality 


cph{x)f{x)  dx 


+ 


(ph  ix)  (fix)  -  f  (x))  dx 


+ 


(ph  ix)  f  (x)  dx  -  f  ix') 

Jr3 


by  the  triangle  inequality 


<  ^  ^  ^  by  inequalities  (4.30),  (4.31),  and  (4.32) 


□ 

Corollary  4.2.  If  the  time  discretized  space  homogeneous  Boltzmann  equation  (3.2)  has 
a  non-negative  bounded  solution  f,  then  the  solution  f  (4.9)  of  the  DMC-KDE  method 
converges  pointwise  to  f. 

Proof.  The  proof  follows  naturally  from  Corollary  4.1.  □ 

The  current  work  represents  the  first  time  these  forms  of  convergence  have  been 
proven  for  a  stochastic  particle  method,  and  illustrates  that  DMC-KDE  transforms  DSMC 
from  a  Boltzmann  simulator  to  a  Boltzmann  solver,  allowing  direct  evaluation  of  the 
velocity  density  function. 
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4.4  Application  to  Bobylev,  Krook,  and  Wu  Problem 

To  provide  a  numerical  demonstration  of  the  DMC-KDE  technique,  the  method  was 
applied  to  the  Bobylev-Krook-Wu  problem.  As  with  Nanbu’s  DSMC  technique,  it  is 
possible  to  compute  the  normalized  moments  by  substituting  the  DMC-KDE  distribution 
function  into  (2.61) 


Zn  ( t ) 


_ 1 _ v_L  f 

(2m  +  1)!!  Zj  1 v/3  J„ 3 


K 


v  ~  Vi  ( t ) 
h 


|2« 


dv 


(4.33) 


from  which  closed  form  expressions  for  the  moments  can  be  derived.  The  first  three 
moments  are  given  by 


z  i  (!) 


Z2  (0 


=  3 ki^  +  lh2 


P  i—  1 


15  A, 


p  i= i 


Z3(')  =  iskZii‘«+ 


p  i=  i 


107 
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One  may  therefore  see  that  the  Nanbu  solution  and  the  DMC-KDE  solution  for  the 
normalized  moments  vary  only  by  a  term  involving  h2'\  Since  h  0  as  Np  — >  oo  these 
terms  become  negligible  as  the  number  of  simulated  particles  become  large  and  are  rapidly 
diminishing  for  higher  order  moments. 

Figure  4. 1  presents  a  comparison  of  the  first  three  normalized  moments  of  the  density 
function  for  the  case  in  which  Np  =  100.  The  results  have  been  averaged  over  an  ensemble 
of  600  runs.  Note  that  like  the  DSMC  results,  there  is  a  significant  amount  of  variance 

present  in  the  solution  and  that  the  solution  tends  to  diverge  more  significantly  for  the  higher 
order  moments.  Additionally,  note  that  the  moments  are  offset  by  the  hr"  term,  which  is 
most  visible  for  the  first  moment. 
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Figure  4.1:  Normalized  Moments  of  Bobylev-Krook-Wu  Solution  (Np  =  100,  600 
run  ensemble) 


Figures  4.2  and  4.3  illustrate  the  same  normalized  moments  for  the  DMC-KDE 
solution  of  the  Bobylev  problem  for  an  ensemble  of  400  samples.  Notice  that  in 
comparison  to  Figures  3.2  and  3.3,  the  moments  are  offset  by  an  amount  proportional 
to  h2n.  As  predicted  theoretically,  the  offset  diminishes  rapidly  with  increasing 


Figure  4.2:  First  Normalized  Moment  with  Varying  Np  (400  run  sample) 
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1.2] 


Figure  4.3:  Second  Normalized  Moment  with  Varying  Np  (400  run  sample) 

moment  orders  and  as  the  number  of  simulated  particles  is  increased.  Strictly  speaking, 
however,  notice  that  neither  method  performs  particularly  well  in  capturing  the  evolution 
of  the  higher  order  moments.  The  advantage  of  DMC-KDE  is  that  it  allows  for  direct 
computation/visualization  of  the  distribution  function,  but  when  only  examining  moments 
this  advantage  is  not  observable. 

To  better  exhibit  this  advantage  consider  the  L1  HR3)  error  between  the  approximate 
and  actual  solution  defined  by 

Error  (f)  =  f  |/(v)  -  /(v)|  dv  (4.34) 

JR3 

As  the  Bobylev-Krook-Wu  solution  is  bounded  and  continuous  for  all  t,  Corollary  4.2 
guarantees  that  /  (v)  — »  /(v)  for  all  v  6  M3.  Since  /  is  also  bounded  and  continuous 
it  is  trivial  to  prove  that  in  this  case  convergence  in  the  L1  (M3)  norm  is  also  guaranteed. 
Figure  4.4  presents  the  results  for  the  L1  (lR3J  error  as  it  varies  with  the  number  of  simulated 
particles  for  a  1000  run  sample. 

Notice  that  the  convergence  is  fairly  monotonic  with  some  exception  at  larger  values 
of  t.  As  stated  in  the  previous  section,  this  form  of  convergence  has  never  previously  been 
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Figure  4.4:  L1  (m3)  Error  for  DMC-KDE  Bobylev-Krook-Wu  Solution  with 
Varying  Np  (  100  run  sample) 


proven  for  a  stochastic  particle  scheme  for  the  Boltzmann  equation.  The  present  work  has 
demonstrated  that  by  applying  principles  from  kernel  density  estimation  such  schemes  can 
be  taken  from  the  realm  of  Boltzmann  simulators  to  the  realm  of  Boltzmann  solvers. 
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V.  A  General  Approach  for  Distributional  Monte  Carlo  (DMC)  Methods  for  the 
Space  Homogeneous  Boltzmann  Equation 

As  observed  in  the  previous  chapter,  the  DMC-KDE  approach  allowed  for  stronger 
forms  of  convergence  than  traditional  DSMC  and  provides  a  means  for  recovering  the 
velocity  density  function.  The  fixed  functional  form  of  the  particle  velocity  density 
functions,  however,  precludes  it  from  being  considered  a  fully  distributional  method  as 
envisioned  by  the  author.  A  Distributional  Monte  Carlo  method  for  the  Boltzmann  equation 
is  a  particle  simulation  which  employs  a  non-singular  representation  of  the  velocity 
density  function  by  allowing  each  simulated  particle  to  posess  a  velocity  density  function 
throughout  the  simulation  instead  of  just  a  single  velocity.  Collision  outcomes  for  a  given 
pair  are  determined  by  computing  an  approximate  space  homogeneous  relaxation  of  the 
velocity  density  functions  of  the  two  simulated  particles  participating  in  a  collision. 

From  a  physically  intuitive  point  of  view,  this  approach  makes  sense.  As  each 
simulated  particle  represents  a  very  large  number  of  actual  particles,  kinetic  theory  tells  us 
that  intermolecular  collisions  occurring  within  the  collection  will  drive  the  density  function 
toward  equilibrium.  What  is  perhaps  somewhat  suprising  is  that  it  can  be  shown  that  such 
an  approach  arises  naturally  from  the  time  discretized,  space  homogeneous  Boltzmann 
equation  when  nonsingular  particle  distributions  with  arbitrary  forms  are  employed. 

In  this  case,  the  velocity  density  function  of  the  gas  is  modelled  by 

1  Np 

/ (v)  =  Tf  y  g»  (v)  (5.1) 

iyp  i=  1 

where  is  the  velocity  density  function  of  the  ith  simulated  particle.  No  particular 
functional  form  of  g,  is  assumed  in  the  general  approach.  Let  f,g  be  density  functions 
such  that  /(v)(l  +  ||v||)  6  Ll  (lR3),  g(v)  (1  +  ||v||)  6  L1  (M3),  /  (v)  log  (/  (v))  6  L1  (M3), 
and  g  (v)  log  (g  (v))  e  L1  (iR3).  Then  by  Theorem  2.3  the  space  homogeneous  Boltzmann 
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equation  has  a  unique  solution  over  the  time  interval  At  for  the  initial  condition  f0  = 
+  g).  Denote  this  solution  by  Q  (/,  g ). 

To  illustrate  the  mathematical  motivation  for  treating  collision  interactions  as  a 
relaxation  problem,  recall  the  weak  transition  kernel  formulation  of  the  space  homogeneous 
Boltzmann  equation  from  Chapter  2 


(f)  (v)  fk+]  (v)  dv=  f  f  KVyW  (cp)  f  (v)  f  (w)  dwdv 

Jr3  Jr3  Jr3 


(5.2) 


By  definition  we  have 


f 

Jr3 


<P  (v)  Q  (/,  g)  dv 


Kv,w  {(p)  (/  (v)  +  g  (v))  (/  (w)  +  g  (w))  dwdv 

Kv,w  {(p)  (/  (v)  f(w)+f  (v)  g  (w)  +  g  (v)  /  (w)  + 
g  (v)  g  (w))  dwdv 

&<j,f)dv  +  2 


HL 
L 


Ky,w  (0)  /  (v)  g  (  w)  dwdv+ 


@{g,g)dv 


(5.3) 


Rearraging,  one  obtains 

Kv,w  (cp)  f  (v)  g  (w)  dwdv 


f 

Jr3 


(p  (v)  ( W  (/,  g)  -  X-Q  (/,  /)  -  (g,  g)  I  dv  (5.4) 


Allowing  fk  to  have  the  form  of  equation  5.1,  upon  direct  substitution  into  Equation  5.2, 
one  obtains 


/ 

Jr3 


cp(v)fk+1  (v)dv  = 


i  N  N  p  p 

1  VVf  f 


KVtW  {(p)  g-  (v)  g)  (w)  dwdv 


(5.5) 


N2  “  ■“  Jr3  Jr3 

Note  that  each  particle  pair  contibutes  a  term  similar  in  form  to  Equation  5.4  to  the  solution, 
therefore  one  obtains 

(p(v)fk+1  (y)dv 


r 

Jr3 


&£££,* M  (2@  (*?•  *5)  -  \6  &  *?)  -  \s  (*!•  4))  * 

(5.6) 
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This  demonstrates  that  the  solution  at  the  next  time  step  can  be  written  as  a  summation 
of  space  homogeneous  relaxation  solutions  for  the  particle  pairs.  This  expression  can  be 
further  simplified  to 


(v)dv  = 


2 

N* 


(p(v)^[g-,gkj)dv+  ^ 


i£^Msu) 


dv 


(5.7) 


In  this  form,  the  evolved  solution  is  clearly  decomposed  into  two  contributions.  The  first 
term  on  the  right  side  represents  space  homogeneous  relaxation  among  the  particle  pairs 
while  the  second  term  represents  the  self-evolution  of  the  individual  particle  distribution 
function  due  to  intermolecular  collisions  only  between  the  actual  particles  it  represents.  In 
the  special  case  where  the  g\  terms  are  Maxwellian,  Q  (gf, g.)  =  gf,  and  thus 


f 

J  R3 


<p(v)fk+l  (v)  dv  = 


N 2 
2 

ac 


§f?l 


N 


dv 


«.=*  j 

'*J  j+i 
N  N 


iU 


<P  (v)  gkj  (v)  dv 


<p(v)f(v)dv  (5.8) 


«.=*  j 
>*J  j*i 


which  is  somewhat  analguous  to  the  form  obtained  by  Nanbu  in  Equation  3.20  in  that  it  is 
decomposed  into  a  binary  interaction  term  and  the  original  density  function. 

In  principle,  any  stochastic  scheme  which  emulates  the  density  function  on  the  right- 
hand  side  of  Equation  5.6  should  be  a  valid  simulation  scheme  for  the  space-homogeneous 
Boltzmann  equation.  Within  the  Distributional  Monte  Carlo  framework,  one  such  choice 
is  as  follows 


•  Choose  a  set  of  velocity  density  functions  for  which  the  probability  measure, 
PN  on  M3  defined  by 

l  N  r 

Pn(A)  =  mTj  JAgi(v)dv  (5'9) 

converges  in  law  to  P0  defined  by 


P o  (A)  =  f  f0  (v)  dv 

J  A 


(5.10) 
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where  /0  is  the  initial  condition  for  the  Boltzmann  equation. 

•  Choose  a  sequence  of  equidistributed  random  numbers  {r,},  r,  6  [0,1]. 

•  Define  the  indices  of  collision  partners 

C  (i,  N)  =  LfVrJ  +  1  (5.11) 

with  velocity  density  functions  given  by 

hi  (v)  =  gC(i,N)  (v)  (5.12) 


Compute  the  set  of  velocity  density  functions  {pi)?=l  where  p,  is  given  by 


Pi  =  2 &(gi,hi)  -  ^ 0(gi,gi )  -  i @(hhhi ) 


The  overall  density  function  is  given  by 


i  N 


1=1 

The  sequence  of  probability  measures  {P^l  on  defined  by 

Pn(A)  =  J'  f(v)dv=^Xjf  Pi  (v) dv  for  any  A  c  K3 

will  be  shown  to  converge  in  law  to  P\  defined  by 

Pi  (A)  =  J  f i  00  dv 

as  iV  — >  oo.  Here  /i  is  the  Boltzmann  solution  at  time  t  =  At 


(5.13) 


(5.14) 


(5.15) 


(5.16) 


Repeat  the  process  with  f\  as  the  new  initial  condition. 
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5.1  Proof  of  Convergence  of  DMC  Approach 

In  developing  this  approach,  it  became  necessary  to  develop  a  generalization  of  one 
of  Babovsky’s  theorems.  Namely,  the  following  theorem  provides  a  sufficient  condition 
for  the  convergence  of  any  numerical  scheme  for  the  Boltzmann  equation;  that  is,  no 
assumptions  are  made  regarding  the  use  of  a  point  measure  approximation  to  the  density 
function. 


Theorem  5.1.  Let  fk  be  the  solution  of  the  time  discretized  space  homogeneous  Boltzmann 
equation  (3.2)  at  t  =  4,  and  let  Pk  be  the  probability  measure  defined  by 


Pk  (A) 


I 


fk  (v)  dv  for  any  A  c 


(5.17) 


Define  P2k  to  be  the  probability  measure  on  D  X  l3  X  R3  given  by 


P2k(B)  =  I  I  f  (v)  f  (w)  dxdwdv  (5.18) 

Id 


If  a  sequence  of  probability  measures  |P2WJ  converges  to  P2k  as  N  co,  then  the  sequence 

of  probability  measures,  { P2N  o  VF  1  j,  on  M3  converges  to  Pk+1  in  law  as  N  — >  oo. 

Proof  Since  'P  is  absolutely  continuous  and  P2N  — >  P2k  in  law  as  N  — >  oo,  we  have  that 
P2n  o  T7”1  — >  P2k  o  'P_1  in  law  as  N  — »  oo  by  Lemma  3.4.  Take  any  f  6  Q,  (M3).  Then, 

j  f(v)d(P2k  o¥~‘)  =  j  (x,v,w))dP2k  by  Lemma  3.3 
JR3  V  7  JR3 

=  r  r  0  (T7  (v,  v,  HO)  /A  (v)  fk  (w)  dxdwdv  by  (5 . 1 8) 

JR3  Jr3  Jd 

nKV'W  ((f))  fk  (v)  fk  (w)  dxdwdv  by  (3.35)  and  (3.36) 

3 

=  f  <P  (v)  fk+1  (v)  dv  by  (3.32) 

Jr3 

Thus,  P2k  o  lF  1  =  Pk+1  and  hence,  P2N  o  VF  1  converges  in  law  to  Pk+l  as  N  — >  oo.  □ 


Thus,  any  scheme  which  can  be  shown  to  be  constructed  in  such  a  way  that  the  above 
theorem  is  satisfied  will  be  convergent. 

The  following  results  will  be  needed  for  the  convergence  proofs  that  follow. 


91 


Theorem  5.2.  (Kolmogorov’s  Convergence  Criterion )  [60]  Suppose  { Xn }  is  a  sequence  of 
independent  random  variables.  If 


YjVar(x})<  oo  (5.19) 

i=i 

then 

OO 

^  (Xj  -  EXj )  converges  almost  surely.  (5.20) 

l=i 


Lemma  5.1.  (Kronecker’s  Lemma)  [60]  Given  two  sequences  { xk }  and  {ak}  such  that 


xk,  ak  6  M  and  0  <  ak  T  °°.  If 

oo 

/  —  converges,  (5.21) 


then 


lim  an  1 

n—>oo 


n 


^xk  =  0 
k=  1 


(5.22) 


While  the  traditional  strong  law  of  large  numbers  typically  is  stated  as  applying  to 
only  sequences  of  independent,  identically  distributed  random  variables,  by  combining 
Theorem  5.2  and  Lemma  5.1  one  can  derive  the  following  Corollary  for  independent 
sequences  of  random  variables  (not  neccessarily  identically  distributed),  which  is 
sometimes  labelled  as  Kolmogorov’s  Strong  Law  of  Large  Numbers. 


Corollary  5.1.  (Kolmogorov’s  Strong  Law  of  Large  Numbers)  [60]  Let  {X„}  be  a  sequence 
of  independent  random  variables  satisfying  E\X„ )  <  oo.  If  { b„ }  is  a  monotonic  sequence 
such  that  bn  T  oo,  and  if 


(5.23) 


then 


ES, 


0  almost  surely  as  n 


(5.24) 


where  Sn  =  Xl'=i  Xk. 
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This  result  allows  one  to  prove  the  following. 


Theorem  5.3.  Let  fk  be  the  solution  to  the  time  discretized  space  homogeneous  Boltzmann 
equation  (3.2),  and  define  P  to  be  the  probability  measure  on  M3  given  by 


P(A)=  I  f  (v)  dv  for  any  A  c 


(5.25) 


Define  the  probability  measure,  P  on  M3  X  M3  by 


P  (B)  =  /  (v)  /  (w)  dvdw  for  any  B  cRjx 


3  v/  TO  3 


(5.26) 


Let  PN  be  the  probability  measure  defined  on  M3  by 


Pn(A )  =  - 


■HI 


gi  (v)  dv  for  any  A  c 


(5.27) 


where  is  a  set  of  density  funtions.  Also  define  the  probability  measure  Pn  on  M3  xl3 


Pn(B) 


(v)  hj  (w)  dvdw  for  any  B  C  M  X 


3  v  m3 


(5.28) 


where  h,  are  defined  as  in  (5.12).  IfPN  converges  to  P  in  law  as  N  — »  oo,  then  PN  converges 
in  law  to  P  as  N  oo. 


Proof.  Define  the  distribution  functions  associated  with  P  and  P  to  be 


F(v)=  f  f(v')dv'  (5.29) 

kJv'<V 

G(v,w)=  f  ^  f  (v')f  (w')dv'dw'  =  F(v)F(w)  (5.30) 

kJv'<V  *Jw'<w 


respectively.  Also  define  the  distribution  functions  associated  with  PN  and  PN  to  be 

T  1  N 

Fn(v)=  —  J^giivfidv' 

Jv'<v  M  i=l 

r  r  i  A 

Gn(v,  w)  =  I  I  —  Y'gi(v,)hi(w,)dv,dw' 

Jv'<v  Jv'<w  W  /=1 


(5.31) 


(5.32) 
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Take  any  ( a,b )  e  M3  x  M3.  We  will  prove  GN(a,b )  — >  G(a,b )  as  iV  — >  oo,  and  thus  TV 
converges  to  P  in  law  as  iV  — >  oo.  Define  the  particle  distribution  functions 

^(v)=  f  gi(y')dv'  (5.33) 

Jv'<v 

Note  that 

1  * 

FN(v)=^YjTi{V)  (534) 

i=l 

Let  Q„  =  (0, 1).  Then  (Q„,  Sn,  A)  is  a  probability  space,  where  Sn  is  the  standard  cr-algebra 
on  (0, 1)  and  A  is  Lebesgue  measure.  Define  Q  to  be  the  cartesian  product 

N 

n  =  n*  (535) 

n=  1 

The  elements  of  Q  are  thus  AMuplcs  {m„}^=1  with  a>„  e  Q„.  Defining  nm  to  be  the  natural 
projection  of  Q  onto  Qm  (i.e.  nm  ( ! j  =  a>m),  let  S  be  the  smallest  cr-  algebra  of 
subsets  of  Q  containing  all  sets  n~l  (A)  for  all  n  and  all  A  c  Q„.  Here,  n~l  (A)  denotes 
the  pre-image  of  A  (namely,  n~]  (A)  :=  (wefl:  n(co)  6  A  c  D„})  .  Then  (D,  S,  Prob)  is 
a  product  probability  space  with  Prob  =  AN.  Define  the  sequence  of  random  variables 
Xi  :  Q  — >  [0, 1]  by 

X,(co)  = 


where 

Note  that  by  Equation  5.32, 


I  I  gi  (v)  hi  (w)  dvdw 

*Jv<a  *Jw<b 

I  gi  (v)  dv  I  gC(i,N)(yv)dw 

*Jv<a  w<b 

T,  ( a )  I  gcu,N)  (w)  dw 

*Jw<b 

73  (fl)  P'c(i,N)  ( b ) 

Ti(a)Yi(aj) 


(5.36) 


Yt  (to)  =  Tam  (b) 


(5.37) 


1  N 

Gn  ( a ,  b)  =  —  V  X, 


iV 


(5.38) 
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Unlike  in  Babovsky’s  proof  for  Nanbu’s  scheme,  the  variables  {X,}  are  independent,  but 
not  identically  distributed.  Therefore,  we  cannot  employ  the  same  arguments  based  on 
the  Central  Limit  Theorem  that  Babovsky  utilized,  but  rather  will  employ  Kolmogorov’s 
Strong  Law  of  Large  Numbers.  To  begin,  recall  £X,  =  T\  ( a )  EYh  thus  consider 


EYi  = 


(5.39) 


Recall  from  (5.11)  that  Tc(i,N)  (b)  can  be  written  piecewise  over  the  intervals  for 

j  =  l,..., N,  namely, 


% ca,N )  (b)  -  Tj  (b)  for  <ol  e 


j  ~  1  j_ 
N  ’  N 


(5.40) 


Substituting  (5.40)  into  (5.39), 

u/v 


-a  IN 


EYi  = 


T\  (b)  dA  +  I  T2(b)dA  + 

0  Jl/N 

hi™ 

i=  1 

EN(b) 


+  /‘ 


TN{b)dA 


{N-D/N 


(5.41) 


Thus, 


EXi  =  Ti{a)F  (b) 


(5.42) 
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Also,  recall  Var(Xj)  =  (T,  (a))'  Var  (Y,).  Thus  we  have 


Var (Yd  =  E(Y'f)-(EYl)2 

=  f  Yf  (oj)  dAN  -  ( F  ( b ))2  by  (5.39) 

Jn 

=  ( W)  (b)f  dAN  -  ( F  (b))2  by  (5.37) 

Jni  J  fl,  J  n^ 

=  f  (rc(i,N)(b))2dA-(FN(b))2 
Jn, 


Jr*l/N  r*A 

C Ti(b)?dA+  ( T2(b))2dA  +  ---  + 

0  Jl/iV 

f  (T/V  (Z>))2  -  (F*  (Z>))2  by  (5.40) 

J(V-l)/iV 
N 

,  (TAhYr  -(F.AhYc 
N 


<L/N 


r=l 

<  \-(FN(b))2 

<  1 


(5.43) 


Consider 


N  j  V 

ZpFar«)  =  z 

A=1  fc=l 


w  on  (a))2 


*=1 

V 


*  X 


k2 

(Tk(a)f 


k=  i 
v 


F 


Var(Yk) 


by  (5.43) 


< 


yl 

hkl 


<  oo 


Therefore,  having  satified  the  requirements  of  Corollary  5.1,  we  conclude 

S  N  ~  ES  N 


N 


0  almost  surely 


where  S N  =  Yj?=i  %i  =  NGn  ( a ,  b).  Thus, 


iV 


ESn  =  ^  70(a)  Fv(F 

i=l 

=  NFN(a)FN(b ) 


(5.44) 
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By  Equation  5.44  we  have  that 


Gn  ( a ,  b )  — >  Fn  ( a )  F#  (&)  almost  surely.  (5.45) 

Recall  that  by  assumption  FN  converges  to  F.  Thus, 

Fn  (a)  Fn  ( b )  ->  F  (a)  F  ( b )  =  G  (a,  b)  (5 .46) 

Combining  this  with  Equation  5.45  yeilds 

GN(a,b)  — >  G(a,b)  almost  surely  (5.47) 

for  all  (d,l))el3xl3,  which  implies  converges  to  P  in  law  as  ./V  — >  oo.  □ 


We  are  now  ready  to  prove  our  main  result. 


Theorem  5.4.  Let  fk  be  the  solution  of  the  time  discretized  space  homogeneous  Boltzmann 
equation  (3.2)  and  define  the  probability  measure  Pk  on  M3  by 


1 


Pk  (A)  =  f  (v)  dv  for  any  A  c 


(5.48) 


Given  a  set  of  density  functions  {gi}^1  define  the  sequence  of  probability  measures  j  P'N  J  on 

M3  by 

i  A  r 

(5.49) 


i  N  r 

pn(a )  =  JASi(v)dV 


Let 


Pi  =  2 QigM  -  -Q(gi,gi)  -  -Q  (hj,  hd 


(5.50) 


with  hj  as  given  in  (5.12).  Define  the  sequence  of  probability  measures  { P'f 1 }  on  R3  by 

Pi  (v)  dv  (5.51) 

If  PkN  converges  to  Pk  in  law  as  N  — »  oo,  then  Pkf 1  converges  to  Pk+l  in  law  as  N  — >  oo 
almost  surely. 
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Proof.  Take  any  (p  e  Cb  (R-).  Then, 


r  c  1  \ '  1 

I  <P(v)=  (j>  (v)  —  V  2Q  [gh  gC(i,N)]  (v)  -  -Q  [gi,  gi 

Jr3  Jr3  [N  ~i  2 


(5.52) 


2  &  [SauN),  gcuv)]  (v)  dv 


Employing  (3.32),  one  may  obtain, 


f  4>(y)Q  [gi,  gC</.,V,J  (y)dv=\  f  f  K,,w  (0)  [gi  (v) gi  (w)  +  gi  (v)  gca.N)  (w)  + 

Jr3  4  Jr3  Jr3 

gi  (w)  gc(/,A0  (v)  +  gca.N)  0)  gc(/,w)  O)]  dvdw 

=  7  f  f  KV'W(4>)gi(v)gi(w)dvdw  + 

4  Jr3  Jr3 

-  I  I  KVtW{f)gi(y)gca,N)(w)dvdw  + 

4  Jr3  Jr3 

-  Kv^w  (0)  gi  (w)  gca,N)  (v)  dvdw  + 

4  Jr3  Jr3 


i  r  r 

4  Jr3  Jr3 


^v,h-  (0)  gca.N)  (v)  gcow)  (w)  dvdw  (5.53) 


Recall  that  (0)  =  ^vr  v  (0),  thus  (5.53)  reduces  to 


f  <P(v)Q  [gi,  gca.N )}  (v)  dv=  \  [  f  KVt w  (cp)  gi  (v)  y,  (w)  dvdw  + 

Jr3  4  Jr3  Jr3 


^v,w  (0)  gi  (v)  gC{i.N)  (w)  JvJw  + 

I3 

Kv,w  W  gC(i,N)  (v)  gca.N)  (w)  JvJw 


Recalling  the  definition  of  Q.  we  have 


f  0  (v)  ^  [gi,  gc(i.V)]  (v)  dv  =  j  f  f  (v)  £  [gi,  gf]  (v)  Jv  + 

Jr3  4  Jr3 

-  I  I  KV'W((p)gi(v)gc(i,N)(w)dvdw  + 

^  Jr3  Jr3 

<P  (v)  £  [gc(i.v),  gca.N)]  (v)  dv 

Jr3 


(5.54) 


Rearranging,  yields 


nKv,w  (0)  gi  (v)  gca.N)  (w)  dvdw  =  2  <p  (v)  £  [g;.  gC(i, ao]  (v)  dv  - 

3  Jr3 

l.  f  f(v)Q  [g„  gi]  (v)  dv  -  ^  f  (p  (v)  Q  [gCa.N),  gc(iW)]  (v)  dv 
4  Jr3  4  Jr3 


(5.55) 


98 


Substituting  (5.55)  into  (5.52)  yields, 

/  =  fiZf  f 

Jr 3  A  JK3  Jr 


Kv,w  (0)  gi  (v)  gc(i,N)  (W)  rfvdw 


■  mu. 


(p  (VF  (jc,  v,  w))  gi  (v)  hi  (w)  dxdvdw  (5.56) 


by  Lemma  3.35.  Denote  by  P2W  the  probability  measure  onDxK3xK3  given  by 


P2n(B) 
Recalling  (5.56)  we  have 


r  i  N 

1*5 


gi  (v)  /?,  (w)  dxdvdw  for  any  Bcflxtfx 


(5.57) 


f  cf>(v)dPkN+l  =  f  f  f  0(W(x,v,w))dP2Nby(5.57) 

Jr3  Jd  Jr3  Jr3 

=  I  (p(v)d(P2N  o  'F-1)  by  Lemma  3.3 
Jr3  v  ’ 

Thus,  we  have  shown  that  P1^  1  =  /J2,.v  o  VF  1 .  Denote  by  P2  the  probability  measure  on 
Dxl3xl3  given  by 


(5.58) 


P2  (. B ) 


I 


f  (v)  /  (w)  dxdwdv  for  any  B  c  D  x  R3  x 


Also,  denote  by  PN  the  probability  measure  on  M3  x  M3  given  by 


Pn(C) 


=  iv  f 

NifJc 


gi  (v)  hi  (w)  c/vc/vv  for  any  C  c  M3  x 


(5.59) 


(5.60) 


Likewise,  denote  by  P  the  probability  measure  on  M3  x  M3  given  by 


Pn(C) 


-X 


/  (v)  /  (w)  c/v’c/vi’  for  any  C  cl3x 


(5.61) 


Consider, 


lim 

N — ^CX) 


m(p  (x,  v,  w)  dP2N  =  lim  f  f  f  (p  (x,  v,  w)  X-  V  gi  (v)  ht  (w)  dxdvdw 

3  JV^oo  JD  JR 3  JR3  A 


=  lim 

N^oo 


'  D 


-  Iff 
~  If  f 

Jd  Jr3  Jr3 

=  ff  f 

Jd  Jr3  Jr3 


<j p  (x,  v,  w)  dxdPp/  by  (5.60) 


<p  (x,  v,  w)  dxdP  a.s.  by  Theorem  5.3 


(p  ( x ,  v,  w)  f  (v)  /  (w)  dxdvdw 


(p  ( x ,  v,  w)  dP2 
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Thus,  P 2n  converges  in  law  to  P2  as  N  — »  oo.  Therefore,  by  Lemma  5.1  we  have  that 
P^+1  =  P2n  o  'F_1  converges  to  Pk+l  in  law  as  iV  — >  oo.  □ 


While  the  previous  theorem  proved  convergence  in  law  to  the  time  discretized  space 
homogeneous  Boltzmann  solution,  the  following  Corollaries  prove  convergence  in  solution 
for  solutions  which  are  L°°  (iR3)  or  bounded. 

Corollary  5.2.  If  the  time  discretized  space  homogeneous  Boltzmann  equation  (3.2)  has  a 
non-negative  solution  f  e  L°°  (iR3),  then  the  solution  f  of  the  DMC  method  (5.14)  converges 
in  L°°  (R3)  to  f  at  each  time  step.  That  is, 


lim  /-/II  =0 

M  J  J  I  loo 


Proof.  Take  any  e  >  0.  Since  /,  /  €  L°°  (lR3),  there  exist  BUB2  e  M+  such  that  | /  (x)|  <  B\, 
and  |  /( x)\  <  B2  almost  everywhere.  Let  S  \  and  S2  be  the  sets  of  zero  measure  over  which 
these  inequalities  do  not  hold  for  /  and  /  respectively.  Let  S  =  S  i  U  S  2  ■  For  any  x'  6  M3  -S , 
define 


<ph  O)  = 


1 


x/2jrh3 


exp 


lx  -  x 


h2 


where  /i  e  1+.  By  Lemma  4.1,  (f,  is  a  delta  family  as  h  —>  0+  centered  at  x' .  Therefore, 
there  exists  H\  >  0  such  that 


/ 

J  R3 


(ph  (x)  /  (x)  dx  -  /  (x  ) 
for  all  h  <  H\ .  Likewise,  there  exists  H2>  0  such  that 

fh  (x)  /  (x)  dx-  f(x) 


r 

Jr3 


6 

<  - 

3 


6 

<  - 

3 


(5.62) 


(5.63) 


Note  that  //,  is  everywhere  continuous  and  |//,  (x)|  <  f°r  x  6  R3-  Thus  by 

Theorem  5.4,  there  exists  N  >  0  such  that 


J"  (ph  (x)  [f(x)  -  f  (x))  dx 


e 

<  - 

3 


(5.64) 
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for  all  Np  >  N.  Choose  H  <  min  {HUH2},  then  for  any  Np  <  N  we  have, 


\\f~f\L  =  esssup|/-/| 

=  sup  \f  (x')  -  f  (x')\ 

x'eR  3-S 


<  sup 


f{x')~  f  (ph(x)f(x)  dx 

X'6K3-S  I  JR3 

f  (ph  (x)  f(x)  dx-  fix') 

J  R3 


+ 


by  the  triangle  inequality 


<  sup 

x'eR3-S 


fix')-  " 

Jr3 


(i h  (x)  f  (x)  dx 


+ 


(ph  (x)  (/  (x)  -  f  (*))  dx 


+ 


<ph  (x)  f  (x)  dx  -fix') 

J  R3 


by  the  triangle  inequality 


<  ^  ^  ^  by  inequalities  (5.62),  (5.63),  and  (5.64) 


=  e 


□ 

Corollary  5.3.  If  the  time  discretized  space  homogeneous  Boltzmann  equation  (3.2)  has  a 
non-negative  bounded  solution  f,  then  the  solution  f  of  the  DMC  method  (5.14)  converges 
pointwise  to  f. 

Proof.  The  proof  follows  naturally  from  Corollary  5.2.  □ 

Note  that  as  mentioned  previously,  the  current  work  represents  the  first  time  such 
forms  of  convergence  have  ever  been  proven  for  a  stochastic  particle  method  for  the 
Boltzmann  equation. 

5.2  Distributional  Monte  Carlo  using  the  Bhatnagar- Gross- Krook  Approximation 
(DMC-BGK) 

Note  that  the  proofs  presented  in  the  prior  section  assumed  the  space  homogenous 
relaxation  problem  between  the  two  particles  was  solved  in  a  way  which  is  consistent  with 
the  Boltzmann  solution.  In  practice,  this  will  necesitate  some  numerical  approximation 
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technique;  however,  so  long  as  the  technique  is  consistent  with  the  Boltzmann  equation 
the  prior  theorems  guarantee  the  DMC  approach  will  converge  in  law  to  the  Boltzmann 
solution. 

A  number  of  appropriate  techniques  could  be  applied  toward  this  end.  Such  techniques 
could  be  either  stochastic  or  deterministic.  Some  possible  choices  include  moment 
methods,  spectral  methods,  linearization  of  the  collision  operator,  or  even  simplified 
stochastic  particle  schemes.  In  practice,  it  would  be  desirable  to  reduce  the  complexity 
of  the  scheme  as  a  large  number  of  simulated  collisions  will  require  computation.  A 
complementing  aspect,  however,  is  that  since  a  large  number  of  collisions  will  be  computed, 
and  each  particle  already  carries  a  significant  amount  of  information  in  its  distribution 
function,  the  actual  collision  calculation  can  be  rather  coarse. 

As  a  demonstration  case,  the  present  work  employs  the  Bhatnagar-Gross-Krook 
(BGK)  equation  [11]  to  compute  the  outcome  of  intermolecular  collisions.  The  space 
homogeneous  BGK  equation  is  a  model  equation  for  the  Boltzmann  equation  which 
replaces  the  complex  collision  integral  with  a  simple  relaxation  form. 

— f(v,t)  =  -v[f(v,t)  -fM(v)\  (5.65) 

ot 

Here,  v  is  the  collision  frequency  and  fM  is  the  Maxwellian  which  possesses  the  same 
energy  and  momentum  as  /.  In  the  case  where  v  is  constant  with  respect  to  molecular 
velocity,  equation  5.65  becomes  a  linear,  ordinary  differential  equation  in  t.  Given  initial 
data  /o,  the  solution  of  equation  5.65  is  given  by 

/  (v,  t )  =  <?“’"./()  (v)  +  ( 1  -  e~vt)  fM  (v)  (5 .66) 

Though  the  BGK  equation  is  not  fully  consistent  with  the  Boltzmann  equation,  various 
BGK  implementations  have  served  as  an  appropriate  approximation  in  many  investigations 
of  rarefied  gases  [33,  78],  including  variance  reduced  DSMC  schemes  [51]  and  therefore 
represents  an  appropriate  choice  to  demonstrate  the  benefits  of  a  distributional  approach. 
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The  above  solution  is  employed  to  approximate  the  solution  Q  of  the  previous  sections. 
Denoting  the  respective  velocity  distributions  as  g  and  h,  the  initial  condition  is  given 
by  /  =  1/2  (g  +  h).  The  Maxwellian  which  corresponds  to  this  density  function  can  be 
determined  by  matching  the  first  and  second  moments. 


Where  the  mean  velocity  u  =  f  vf  (v)  dv,  a\ 


fM(v)  =  a]e-a^-u)2 

8(M;HM|2)13/2 


,  O'  2 


(5.67) 

8  (m2  -  ||w||2),  and 


M,  =  /l|v||2  f(v)dv.  The  space  homogeneous  evolution  of  the  combined  distribution 


function  may  then  be  obtained  utilizing  (5.66). 


-vA  t 


Q  [<?,  h]  (v) 


[g  (v)  +  h  (v)]  +  (l  -  e  rAr)  axe 


QT2  (V-U)1 


(5.68) 


Equation  5.68  is  a  closed  form  expression  for  the  evolution  of  the  joint  distribution 
function  of  the  two  simulated  particles  involved  in  a  collision  interaction.  Although  an 
exact  expression,  one  must  determine  an  appropriate  way  to  represent  the  distribution 
functions  g,h  numerically.  Several  approaches  exist.  For  simplicity,  the  current  study 
represents  the  distribution  function  as  a  set  of  values  evaluated  on  a  three  dimensional  grid 
in  velocity  space.  Although  conceptually  simple,  this  approach  would  be  suboptimal  for 
implementation  in  a  production  code  as  a  large  amount  of  data  is  tied  to  each  particle. 
In  a  parallel  implementation  the  latency  associated  with  convecting  a  particle  from  one 
processing  domain  to  another  could  become  nontrivial  because  of  the  amount  of  data  to 
be  transferred  to  the  new  processing  unit.  A  better  approach  would  be  to  approximate 
the  distribution  function  in  closed  form.  For  example,  if  one  expanded  the  distribution 
functions  utilizing  Hermite  polynomials  (as  suggested  by  Grad  [41]),  Equation  5.68  could 
then  be  used  to  evolve  each  of  the  Hermite  coefficients  in  time.  The  polynomial  coefficients 
are  then  the  only  data  carried  with  each  particle  and  the  latency  associated  with  a  parallel 
implementation  would  be  significantly  reduced. 

Computation  of  the  distribution  function  over  a  discrete  grid  as  described  above 
presents  the  method  with  another  tunable  parameter,  namely,  the  grid  spacing.  Grid 
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refinement  in  velocity  space  is  a  problem  that  is  common  with  other  numerical  techniques 
for  the  Boltzmann  equation  (e.g.  direct  solvers)  but  is  not  currently  well  understood.  For 
this  reason,  the  results  and  techniques  described  in  this  work  are  based  only  on  uniform 
grid  spacing  in  velocity  space. 

Before  continuing,  it  should  be  noted  that  the  method  described  here  represents  a 
hybrid  stochastic-deterministic  scheme.  That  is  to  say,  collision  pair  selections  continue  to 
be  performed  stochastically  using  the  Monte  Carlo  scheme  previosly  discussed.  Collision 
outcomes,  however,  are  computed  deterministically  using  equation  5.68.  This  hybridization 
is  made  possible  by  the  fact  that  particles  possess  distributed  velocities.  If  the  scheme 
restriced  them  to  a  finite  set  of  velocities,  a  stochasic  approach  would  be  required  to  sample 
from  equation  5.68  (such  a  scheme  was  presented  by  the  author  in  Reference  [67]).  The 
hybridization  of  the  current  scheme,  however,  will  be  observed  in  the  subsequent  section 
to  result  in  a  drastic  reduction  in  variance. 

5.3  Application  to  the  Bobylev-Krook-Wu  Problem 

To  provide  a  numerical  demostration  of  the  DMC-BGK  technique,  the  method  was 
applied  to  the  Bobylev-Krook-Wu  problem.  As  the  solution  is  computed  over  a  discrete 
grid  in  velocity  space,  simple  quadrature  was  utilied  to  compute  the  normalized  moments 
as  defined  by  Equation  2.61.  Figure  5.1  presents  a  comparison  of  the  first  three  normalized 
moments  of  the  density  function  for  the  case  in  which  Np  =  100.  The  cell  width  on  the 
velocity  grid  was  chosen  to  be  Av  =  2/3  over  the  region  [-5, 5]  x  [-5, 5]  x  [-5, 5].  Outside 
of  this  region  the  value  of  the  density  function  was  assumed  to  be  zero.  The  results  have 
been  averaged  over  an  ensemble  of  600  runs. 

Note  that  unlike  the  DSMC  results  or  the  DMC-KDE  results,  the  current  method  does 
not  display  large  amounts  of  variance.  This  effect  is  due  to  the  fact  that  collision  outcomes 
are  determined  deterministically  which  was  facilitated  by  the  distributed  particle  velocities. 
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1.1 


Figure  5.1:  Normalized  Moments  of  Bobylev-Krook-Wu  Solution  (Np  =  100, 
Av  =  2/3,  600  run  ensemble) 


Furthermore,  note  that  not  only  does  the  current  solution  display  reduced  variance,  but  is 
also  much  more  accurate  than  the  prior  methods  discused.  As  each  particle  possesses  an 
entire  distribution  function  of  unrestricted  form,  a  tremendous  amount  of  information  is 
represented  with  each  particle  which  allows  for  the  overall  density  function  to  be  calculated 
much  more  accurately  for  the  same  number  of  particles.  Finally,  also  note  that  as  a  result 
of  the  variance  reduction  previously  discussed,  the  scheme  conserves  energy  much  better 
that  the  prior  schemes  as  evidenced  by  the  nearly  constant  first  normalized  even  moment 
Z\. 

Figures  5.2  and  5.3  illustrate  the  results  for  the  first  and  second  normalized  moments 
obtained  by  varying  Np  while  holding  the  cell  size  constant  at  Av  =  2/3  and  using  an 
ensemble  of  400  runs. 

Note  that  in  this  case,  Figure  5.2  shows  that  the  total  fluctuation  in  energy  throughout 
the  simulation  is  less  that  one-tenth  of  one  percent.  In  contrast,  the  results  for  Nanbu’s 
method  (Figure  3.1)  displayed  about  an  order  of  magnitude  greater  fluctuation  in  this 
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Figure  5.2:  First  Normalized  Moment  with  Varying  Np  (Av  =  2/3,  400  run  sample) 


Figure  5.3:  Second  Normalized  Moment  with  Varying  Np  (Av  =  2/3,  400  run 
sample) 


moment.  The  results  for  the  second  moment  in  Figure  5.3  show  very  little  variance  as 
well  and  nearly  monotonic  convergence. 

As  the  entire  distribution  function  is  computed,  like  DMC-KDE  it  is  possible  to 
compute  the  L1  (m3)  error  via  Equation  4.34.  Figure  5.4  shows  how  the  L1  (R3)  error 
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varies  as  a  function  of  the  number  of  simulated  particles  with  Av  =  2/3.  Note  that 


Figure  5.4:  L1  (R3j  Error  for  DMC-BGK  Bobylev-Krook-Wu  Solution  with 
Varying  Np  ( 100  run  sample) 

the  convergence  is  monatonic,  and  the  maximum  error  throughout  the  simulation  is  about 
one  order  of  magnitude  smaller  than  the  DMC-KDE  approach  (Figure  4.4).  Also,  unlike 
the  DMC-KDE  approach,  note  that  the  error  is  nearly  zero  at  t  =  0.  This  is  due  to  the 
fact  that  the  initial  condition  can  be  much  better  represented  with  non-gaussian  particle 
distributions. 

Figure  5.5  presents  the  results  for  the  Ll  (M3)  error  as  it  varies  with  velocity  grid  cell 
size.  The  simulations  employed  50  simulated  particles  and  were  ensemble  averaged  over 
100  runs.  Here,  the  expected  general  trend  is  observed,  namely  that  reduced  cell  size 
results  in  reduced  error,  however  the  results  are  not  necessarily  monatonic.  This  is  likely 
owed  to  the  fact  that  as  domain  was  finite  and  held  constant,  at  some  point  the  error  in 
approximating  the  quadrature  becomes  smaller  than  the  stochastic  variance  introduced  by 
collision  pair  selection  over  a  collection  of  only  50  particles. 
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0.12 


Figure  5.5:  L1  (R3)  Error  for  DMC-BGK  Bobylev-Krook-Wu  Solution  with 
Varying  Av  (100  run  sample) 

The  total  variation  for  the  first  three  moments  as  a  function  of  the  size  of  the  ensemble 
average  is  given  in  Figure  5.6.  Here  Np  =  100  and  Av  =  2/3.  Notice  the  drastic 
improvement  over  Figure  3.4.  The  total  variation  is  very  nearly  equal  to  the  analytical 
value  even  for  a  very  modest  number  of  samples.  This  is  owed  in  large  part  to  the  fact 
that  the  stochastic  step  predicting  collision  outcomes  has  been  replaced  by  a  deterministic 
calculation,  which  due  to  the  ability  of  particles  to  posess  distributed  velocities  accounts 
for  all  possible  collision  outcomes.  This  represents  the  first  use  of  deterministic  collision 
modeling  within  a  stochastic  particle  method. 

It  should  be  noted  that  although  the  numerical  solutions  presented  in  this  section 
compare  well  to  the  Boltzmann  solution,  it  should  not  be  assumed  that  the  scheme 
employing  the  BGK  approximation  is  convergent  to  the  Boltzmann  solution.  Although 
in  many  cases  the  BGK  model  can  be  a  useful  surrogate,  it  cannot  be  directly  derived 
by  simplification  of  the  Boltzmann  collision  integral.  Nevertheless,  the  results  presented 
here  suggest  that  so  long  as  an  appropriate  number  of  collisions  are  calculated,  some 
simplification  of  the  collision  operator  is  tolerable. 
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Figure  5.6:  Total  Variation  of  DMC-BGK  Method  as  a  Function  of  Sample  Size  for 
Np  =  100,  Av  =  2/3 
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VI.  Conclusion 


The  work  presented  herein  has  contributed  a  number  of  important  discoveries  and 
developments  to  the  stochastic  solution  of  the  Boltzmann  equation.  First  and  foremost, 
this  research  has  introduced  the  Distributional  Monte  Carlo  concept  for  the  Boltzmann 
equation.  The  framework  derived  by  the  author  represents  the  first  stochastic  particle 
based  method  which  does  not  employ  a  point  measure  based  approximation  of  the  velocity 
density  function,  but  rather  allows  each  simulated  particle  to  possess  an  entire  velocity 
density  function  of  its  own.  Allowing  the  simulated  particle  velocities  to  be  distributed 
enables  the  modeling  of  collision  events  between  simulated  particles  to  be  treated  as  a  space 
homogeneous  relaxation  problem  over  the  collection  of  actual  particles  they  represent. 
While  physically  intuitive,  the  author  was  able  to  prove  that  such  a  representation  results 
naturally  from  the  weak  time  discretized  space  homogeneous  equation  when  particle 
velocities  are  distributed.  In  addition  to  deriving  the  general  DMC  approach,  it  was 
also  proven  that  the  technique  was  convergent  to  the  solution  of  the  space  homogeneous 
Boltzmann  equation  in  law,  and  in  solution  for  L°°  HR3)  solutions. 

The  distributional  framework  derived  by  the  author  is  quite  general,  and  the  present 
work  established  a  sufficiency  criterion  by  which  any  distributional  scheme  will  converge 
so  long  as  collision  modeling  is  consistent  with  the  space  homogeneous  solution.  This 
facilitates  possible  employment  of  both  stochastic  and  deterministic  schemes.  The  author 
presented  one  such  scheme  based  on  the  BGK  simplification  of  the  collision  operator  which 
employed  deterministic  modeling  of  collision  effects.  Note  that  while  collision  outcomes 
were  determined  deterministically,  collision  pair  selection  was  performed  stochastically.  In 
contrast  to  DSMC  which  performs  both  steps  stochastically,  the  technique  is  observed  to 
have  drastically  reduced  variance.  This  is  due  to  three  sources.  First,  replacement  of  the 
stochastic  process  for  selecting  collision  outcomes  with  a  deterministic  one  tends  to  reduce 
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variability.  Second,  the  use  of  distributed  velocities  and  particle  distribution  functions  of 
general  form  allows  for  the  calculation  of  not  just  one  possible  collision  outcome,  but 
computes  the  effects  of  all  possible  collision  interactions.  Finally,  the  mere  fact  that  particle 
velocities  are  distributed  means  that  each  simulated  particle  carries  with  it  significantly 
more  information  than  it  would  in  a  traditional  DSMC  approach,  allowing  for  the  overall 
distribution  function  to  be  represented  much  more  accurately.  Reduced  variance  in  solution 
translates  to  a  reduction  in  the  number  of  particles  required  to  accurately  simulate  a  given 
problem  which  can  have  a  drastic  impact  on  computational  time  to  solution. 

The  author  also  applied  kernel  density  estimation  to  DSMC  to  develop  a  new  approach 
called  DMC-KDE.  It  was  shown  that  through  the  application  of  kernel  density  estimation, 
the  DSMC  technique  could  be  made  to  converge  not  just  in  law,  but  in  solution  for  L°°  (if3 ) 
and  bounded  solutions  of  the  space  homogeneous  equation.  This  represents  the  first  time 
that  such  forms  of  convergence  have  been  proven  for  a  stochastic  particle  method  for  the 
Boltzmann  equation.  Simply  put,  the  DMC-KDE  approach  turns  a  stochastic  Boltzmann 
simulator  into  a  stochastic  Boltzmann  solver.  This  facilitates  the  direct  computation  and/or 
visualization  of  the  distribution  function,  something  which  is  not  directly  available  to 
Boltzmann  simulators  which  converge  only  in  law,  such  as  DSMC. 

Although  the  particle  velocities  were  allowed  to  be  distributed,  they  were  constrained 
to  be  distributed  by  a  prescribed  form  of  a  density  function.  In  this  case,  particles  were 
allowed  to  possess  velocities  distributed  according  to  a  Maxwellian  distribution.  This  leads 
to  the  physical  interpretation  that  the  actual  particles  represented  by  a  simulated  particle 
do  not  all  possess  the  same  velocity  vectors,  but  rather  are  in  translational  equilibrium 
with  one  another.  Collision  interactions  were  computed  stochastically  and  had  the  effect  of 
shifting  the  mean  of  the  particle  Maxwellians.  The  rules  governing  the  selection  of  the  new 
mean  were  shown  to  be  identical  to  the  rules  Nanbu  used  to  assign  post  collision  velocities. 
Unfortunately,  the  restrictions  placed  on  the  form  and  evolution  of  the  particle  velocity 
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density  functions  eliminate  any  potential  variance  reduction  by  distributing  velocities, 
and  for  the  Bobylev-Krook-Wu  solution  it  was  observed  that  the  total  variation  in  the 
normalized  moments  was  equal  to  that  of  Nanbu’s  scheme. 

The  current  work  has  opened  up  a  number  of  avenues  for  continued  research  in 
Distributional  Monte  Carlo  methods  for  the  Boltzmann  equation.  The  next  logical  step 
for  such  methods  is  generalization  to  the  full  Boltzmann  equation.  The  approach  for  such 
a  technique  would  largely  mirror  DSMC,  that  is  a  grid  will  be  employed  in  physical  space 
to  ensure  collisions  occur  between  simulated  particles  which  are  near  neighbors.  Particle 
advection  and  collision  interactions  will  be  decoupled  via  the  uncoupling  principle,  and 
within  a  cell  the  solution  will  be  presumed  to  be  space  homogeneous.  The  approach 
for  handling  particle  advection  is  not  necessarily  determined,  as  simulated  particles 
now  possess  not  just  a  single  velocity  vector  but  an  entire  distribution  function.  Early 
experimentation  by  the  author  suggests  one  appropriate  choice  is  to  sample  an  advection 
velocity  at  each  timestep  with  which  to  convect  each  particle  from  its  density  function. 

There  is  great  potential  for  distributional  approaches  in  the  nonhomogeneous  case 
which  arises  from  the  fact  that  the  distribution  function  is  directly  computed  and  available 
throughout  the  computation.  When  the  distribution  function  of  a  given  cell  is  Maxwellian 
or  near-Maxwellian  collisional  effects  can  be  neglected.  To  illustrate  this,  consider  the 
problem  of  a  one-dimensional  shock  wave.  The  typical  approach  for  modeling  such  a  case 
is  to  initialize  the  upstream  half  of  the  domain  to  an  equilibrium  distribution  based  on  the 
freestream  conditions.  The  downstream  half  of  the  domain  is  initialized  to  an  equilibrium 
distribution  based  on  the  downstream  conditions  which  can  be  computed  via  the  normal 
shock  relations.  In  the  case  of  a  distributional  approach,  at  t  =  0,  the  particles  in  upstream 
domain  share  a  common  Maxwellian  density  function,  as  do  the  particles  in  the  downstream 
domain.  As  particles  begin  to  advect,  the  cells  near  the  middle  of  the  domain  become  the 
first  to  possess  non-Maxwellian  distributions,  and  therefore  become  the  only  location  where 
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collisions  are  required  to  be  calculated.  That  is  to  say,  collision  calculations  are  clustered 
in  the  area  of  the  shock  wave  itself.  As  the  simulation  continues,  the  region  grows  slightly 
but  the  majority  of  collisions  are  computed  in  the  shock  region.  The  method  is  therefore 
self  adjusting  to  regions  of  non-equilibrium,  in  that  it  requires  little  to  no  computations 
in  equilibrium  regions  and  focuses  its  computational  power  in  regions  which  are  highly 
non-equilibrium. 

In  addition  to  generalization  to  the  full  Boltzmann  equation,  there  are  a  number  of 
potential  research  areas  that  can  still  be  addressed  in  the  space  homogeneous  setting.  First, 
for  the  numerical  results  presented  here,  the  current  effort  employed  a  grid  in  velocity 
space.  This  is  suboptimal  for  use  in  a  production  code  as  each  simulated  particle  carries 
with  it  an  enormous  amount  of  information  which  could  introduce  significant  latency  in  a 
parallel  environment  as  particles  convect  between  processing  domains.  A  better  approach 
is  likely  to  expand  the  particle  velocity  density  functions  using  Hermite  polynomials,  in 
which  case  a  handful  of  coefficients  would  be  the  only  data  required  to  fully  describe  each 
particles  velocity  density  function.  Second,  in  addition  to  the  numerical  representation 
of  the  density  function,  the  exploration  of  methods  to  model  the  solution  Q  used  in  the 
collision  simulation  of  the  Distributional  Monte  Carlo  approach  is  an  area  which  holds  a 
number  of  interesting  problems.  The  framework  proven  by  the  author  made  little  restriction 
on  the  form  of  such  modeling,  allowing  for  both  stochastic  and  deterministic  techniques. 
Although  the  present  work  has  demonstrated  significant  variance  reductions  are  achievable 
with  a  deterministic  approach,  there  is  value  in  pursuing  stochastic  approaches  as  well.  One 
such  approach  is  to  apply  kernel  density  estimation  to  the  particle  distribution  functions 
themselves  (versus  the  overall  distribution).  The  author  has  presented  such  a  technique  in 
the  past  [67],  however,  some  additional  work  is  required  to  make  the  approach  consistent 
with  the  DMC  framework  presented  in  this  effort. 
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While  the  DMC-KDE  scheme  presented  did  not  exhibit  variance  reduction  over 
DSMC,  it  did  establish  an  interesting  approach  for  turning  Boltzmann  simulators  into 
Boltzmann  solvers.  There  are  a  number  of  modifications  that  could  be  explored  which 
could  enhance  the  performance  of  such  an  approach.  As  mentioned  above,  the  author  has 
presented  a  scheme  which  applied  kernel  density  estimation  at  the  particle  density  function 
level  [67]  which  allowed  for  multiple  velocity  samples  per  particle.  This  has  the  effect 
of  unconstraining  the  particle  density  functions  from  a  fixed  form.  Another,  less  drastic, 
approach  might  be  to  allow  the  kernel  bandwidth  to  vary  for  each  particle.  In  this  case, 
particle  densities  would  continue  to  be  Maxwellian  but  each  particle  would  have  a  unique 
mean  and  standard  deviation.  This  is  equivalent  to  allowing  each  particle  to  possess  a 
translational  temperature.  Relations  governing  the  collision  modeling  would  need  to  be 
developed  which  not  only  determined  a  mean  post-collision  velocity,  but  post-collision 
temperature  for  each  simulated  particle. 

When  viewed  as  a  whole,  the  current  effort  has  provided  a  number  of  firsts  in  the  area 
of  stochastic  particle  schemes  for  the  Boltzmann  equation  and  has  opened  a  new  area  of 
research  with  many  interesting  mathematical  challenges  for  future  researchers.  In  addition, 
the  importance  of  hypersonic  and  rarefied  flows  has  ensured  that  techniques  for  modeling 
the  Boltzmann  equation  will  remain  relevant  to  Air  Force  problems  for  the  foreseeable 
future. 
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